dataesr/scientific-paragraphs-categorization huggingface dataset:
To select genes consistently altered in AD brains, considering the several sources of measurement variations in RNAseq experiments  (Van den Berge et al. , 2017) , we decided to focus only on DEGs replicated in at least two independent datasets obtained from related brain areas. This resulted in a set of 734 DEGs detected in both Mayo and MSBB TL (temporal lobe intersection -TLI), and 31 DEGs shared between ROSMAP and MSBB FL (frontal lobe intersection -FLI) (Supplementary table 2). Among TLI DEGs we observed ABCA1 and 2 (ATP Binding Cassette Subfamily A Member 1 and 2), primarily involved in the maintenance of normal brain homeostasis and associated with AD and other neurological diseases  (Abuznait & Kaddoumi, 2012) ; Complement C1R and C1S, involved in the immune/inflammatory response and previously shown to be upregulated in the brain of a 3 × Tg mouse model of AD when Aβ plaques start to accumulate  (Benoit et al.
Again, we used only genes commonly altered in two datasets (TLI or FLI) to avoid inaccurate results associated with the use of large gene sets in functional analysis  (Subramanian et al. , 2005). We found that TLI DEGs were significantly enriched for terms (GO:BP, GO:CC and KEGG) associated with generic biological processes, such as cell-signaling pathways and cell-cell signaling, whereas the small number of DEGs in the FLI were not significantly enriched for any term (Figure  2C ; Supplementary table 3). The limited number of significant gene set enrichment observed in our analysis after inputting DEGs is in disagreement with results reported by  Canchi et al. (2019). This discrepancy can likely be explained by the use of stringent criteria to detect TLI DEGs in our study (only genes detected in at least two independent datasets with FC>1.
Differential transcript usage analysis reveals novel genes associated with AD pathology Gene-level expression analysis lacks the sensitivity to detect possible changes at the transcript-level caused, for example, by alterations in alternative splicing  (Vitting-Seerup & Sandelin, 2017; Yi et al. , 2018). To overcome this limitation, we used differential transcript usage (DTU) analysis to identify additional alterations of gene expression in the AD brains compared to controls. We observed 2,509 and 1,843 genes with differential transcript usage (gDTU) in the temporal lobe of AD brains studied in the Mayo and MSBB datasets, respectively (Figure  3A -B; Supplementary table  1 ). Similar to what we observed for DEGs, a much smaller number of gDTUs were detected in the frontal lobe, both in ROSMAP and MSBB studies (59 and 855 genes with transcripts altered, respectively).
, 2018). Alternative splicing events may have diverse functional consequences for protein expression, such as shifting the frequency of transcripts containing introns (non-coding) or mRNA stability (nonsense mediated decay) or leading to gain/loss of protein domains, intrinsically disordered regions or signaling peptides  (Vitting-Seerup & Sandelin, 2019). Quantification of these consequences revealed some interesting differences between Mayo and MSBB BM36 (Figure  5C ), the two datasets with largest numbers of gDTUs. Whereas in the Mayo dataset, a high number of isoforms showed loss of coding potential and protein domains, in the MSBB BM36 isoforms showed an even distribution of loss and gain of coding potential or protein domains (Figura 5B). These differences could be at least partly explained by the larger number of gDTUs detected in the Mayo compared to MSBB TL (Figure  3 ) and are likely related to the different median read depth of these datasets (Mayo -12.
Differential expression of genes involved in alternative splicing correlates with isoform switches during disease progression To evaluate whether the emergence of gDTUs could be correlated with AD pathology hallmarks, we quantified the total of gDTUs observed at different disease stages in the MSBB dataset using the Braak classification (Figure  6 ; Supplementary Table  5 ). For this purpose, we subdivided samples in three groups: low Braak (0, 1 and 2) -196 samples (clinical diagnosis: 15 AD and 181 controls); mid Braak (3 and 4) -133 samples (clinical diagnosis: 58 AD and 75 controls); and high Braak (5 and 6) -308 samples (clinical diagnosis: 305 AD and 3 controls). Next, we evaluated the number of gDTUs when comparing individuals at these different stages (Figure  6 ). We observed that most gDTUs were detected only while comparing high with either low or mid Braak stages (Figure  6A-D ).
A similar correlation between pathological burden and differential expression of the same 441 splice-related genes was observed in the comparison between all AD versus control subjects of the different datasets (Figure  6F ). Changes in the expression of those genes were hardly observed in FL (low number of gDTUs -Figure  3 ), but were frequent in TL samples (high number of gDTUs -Figure  3 ), albeit to a lesser extent than that observed in the comparison between different Braak stages (likely due to the effects of combining low, mid and high Braak stages in the AD group). Remarkably, the majority of the splicingrelated genes with altered expression in the Mayo dataset was not reproduced in the MSBB BM36 dataset, and vice-versa (Figure  6F ). This could help to explain the dissimilar consequences of alternative splicing events observed in those datasets (Figure  5C ) and suggest that a myriad of proteins could be involved in altered splicing in the AD brains.
Subsets of these genes were exclusively expressed either in microglial cells (14 out of 116), neurons (12), astrocytes (2), oligodendrocytes (6) or endothelial cells (6), suggesting cell-type specific roles for these AD risk/causal genes. Next, we set out to evaluate the differential expression or transcript usage for these genes. Out of the 116 AD risk/causal genes expressed by brain cell types (Figure  8A ), we observed that 54 were also DEGs/gDTUs in at least one of the bulk RNAseq datasets analyzed. Among those genes, 2 were exclusively identified in the FL (Figure  8B ). We, therefore, decided to focus on the 52 AD risk/causal genes identified in the temporal lobe for further analyses. In this region, we identified 27 and 17 DEGs/gDTUs in the MSBB_TL and Mayo datasets, respectively, including some well-characterized AD risk genes, such as ADAM10 (ADAM Metallopeptidase Domain 10), BIN1, CLU (Clusterin) and TREM2 (Triggering Receptor Expressed On Myeloid Cells 2), and the causal AD genes APP, PSEN1 and 2 (Presenilin 1 and 2) (Figure  8A-B ).
Eight genes. were altered in both datasets (Figure  6B , yellow box; 15,38% of overlap) and were selected for further analysis of isoform switch. Using ISAR to identify the isoforms altered in the AD brains compared to controls, we observed some patterns of isoform switch that could have important functional relevance (Figure  8C  and D ). For instance, while BIN1 transcripts ENST00000316724. 9 (NP_647593. 1 -isoform 1) and ENST00000409400. 1 (NP_647600. 1 -isoform 9) were downregulated, transcripts ENST00000393040. 7 (NP_647598. 1 -isoform 6) and ENST00000462958. 5, ENST0000046611. 5 and ENST00000484253. 1 (intron retention) were upregulated (Figure  8C ). This pattern could lead to a decrease of the neuronal specific BIN1 isoform 1 expression  (Zhou et al. , 2014) , given that retained introns are non-coding sequences.
Using western blotting analysis, we confirmed this decrease of BIN1 isoform 1 protein in the frontal cortex and hippocampus of AD brain samples compared to controls (Supplementary Figure  5 ). We also observed isoform switches in the AD causal gene APP with possible functional consequences in neuronal cells. While two APP isoforms were downregulated (ENST00000348990 and ENST00000354192), the isoforms ENST00000346798 and ENST00000357903 were upregulated in Mayo and MSBB datasets (Figure  8D ). Noteworthy, significantly downregulated APP isoforms lack exon 7, which contains the Kunitz protease inhibitor (KPI) domain. KPI is one of the main serine protease inhibitors and increased KPI(+)APP mRNA and protein expression levels have been described in AD brains and are associated with increased amyloid beta deposition  (Tanzi et al.
Gene isoform switches are mostly evident at late stages of the pathology and correlate with altered expression of genes encoding for splicing-related proteins. Using an indirect approach to assign genes to unique cell types, we are also able to map DEGs/gDTUs to unique cell populations of the adult brain, and our results are comparable to previously published scRNAseq data  (Grubman et al, 2019; Mathys et al, 2019). Finally, we show that a subset of AD causal/risk factors such as APP or BIN1 are differentially expressed in the AD brain. Altogether, our work provides a comprehensive description of regional and cell-type specific gene expression changes in the AD brain and suggests that alternative splicing could be an important mechanism for pathological progression. Despite the availability of RNAseq datasets generated from healthy subjects and AD patients  (Allen et al, 2016; Wang et al, 2018; De Jager et al, 2018) , a systematic evaluation of the gene expression changes in the AD brain, as well as comparisons of.
, 2019). They can also help to explain the low number of DEGs identified in scRNAseq data obtained from the frontal lobe  (Mathys et al. , 2019)  compared to a similar study in the entorhinal cortex  (Grubman et al. , 2019). In order to minimize the variability in RNAseq experiments (  Van den Berge et al. , 2017)  we here focused on DEGs (genes with FC>1. 3 and FDR<0. 01 in AD versus control) detected independently in at least two datasets containing samples of similar brain regions (TLI or FLI). These stringent criteria limited the number of DEGs used in subsequent analyses, but still allowed the uncovering of several genes previously associated with AD pathology, such as ABCA1, ABCA2, CALB1, C1R, C1S, GAD1/2, PVALB, REST, SLC32A1, SST, VGF and VIP  (Abuznait & Kaddoumi, 2012; Beckmann et al.
, 2020; Benoit et al. , 2013; Prévot & Sibille, 2020; Zullo et al. , 2019). The reduced number of DEGs in FLI and TLI likely explains our failure to detect functional annotations associated with synaptic transmission and immune response in GSEA, as previously reported  (Canchi et al. , 2019). However, this study analyzed only the ROSMAP dataset and considered genes with FDR<0. 05 as significant, regardless of the fold change, identifying 1,722 DEGs in AD versus control brains. Besides the questionable meaning of DEGs with very small fold changes, the use of such a large set of genes for GSEA can artificially increase the number of significantly enriched functional annotations and is not advised  (Subramanian et al. , 2005). Nevertheless, our failure to detect key functional annotations associated with AD pathology while inputting TLI DEGs is puzzling and could suggest that DEG analysis fails to detect relevant alterations in gene expression in the AD brain.
, 2012; Norstrom et al. , 2010; Yu et al. , 2016). Moreover, gDTUs identified in two independent datasets (TLI), alone or in combination with TLI DEGs, were enriched for key biological processes involved in AD pathogenesis, such as synaptic communication, immune response, inflammation, endocytosis and cell-signaling  (Canter et al. , 2016). Similar gene set enrichment has been described using the analysis of co-expression modules in bulk RNAseq  (Milind et al. , 2020; Wan et al. , 2020)  or DEG analysis DEGs in unique cell types in scRNAseq  (Mathys et al. , 2018; Grubman et al. , 2018). This could suggest that the combination of DEG and DTU to analyze bulk RNAseq is comparable to scRNAseq regarding the sensitivity to detect gene expression alteration in AD brains. In agreement with this possibility, we were able to assign DEGs and gDTUs to unique cell types and confirm the similarities among cell type-specific functional annotations observed in our work compared to previous scRNAseq studies  (Mathys et al.
, 2018)  and suggests that isoform switches in the gene BIN1 are involved in the reduced expression of the main neuronal BIN1 isoform 1 in AD brains. Metadata obtained from each study was used to classify patients into Control and Alzheimer's disease groups (Supplementary table  12 ). Briefly, for the MSBB dataset, we used patients with Neuropathology Category (NP. 1) labeled as "Control" and "definitive Alzheimer". For the Mayo dataset, we used the "Diagnosis" column of the metadata, selecting only "AD" and "Control" patients. For the ROSMAP dataset, we also used the column "Diagnosis" of the metadata, selecting only "Control" (value = 1) and "Alzheimer with no other conditions" (value = 4). In all those datasets, subjects marked as "AD" showed Braak stage values higher than 4. In the MSBB dataset, CDR scores of AD patients were consistently higher than 2.
In the Mayo and ROSMAP datasets, all AD patients had also a definitive diagnosis according to NINCDS criteria. Co-variates such as "Post-mortem interval (PMI)", "RNA integrity number (RIN)", "Age of death" and "Sex" were balanced among the different groups (Table  1 ; Chi-square test, p>0. 05). We used RIN and PMI as covariates in our model to control for possible "batch effects" (linear regression). For detailed information of all individual samples used in this study, please refer to Supplementary table 13.
Differential gene expression analyses Differentially expressed genes (DEGs) were identified using differential gene expression at transcript-level using DESeq2 R library  (Love et al. , 2014 (Love et al. , , 2019)). To facilitate kallisto output import, transcript-level estimated counts, length and abundance was extracted using tximport function  (Soneson et al. , 2016). As described by Michael Love group, transcript-level differential gene expression enhances analysis resolution  (Love et al. , 2019). Using DESeqDataSetFromTximport, a DESeq2 object was created and filtered using rows with sum of all counts bigger than 10. Next, DESeq function was used with default values. Using the results function, we selected all genes with False Discovery Rate (FDR) < 0. 01 and Fold Change (FC) > 1. 3.
We also used RIN and PMI as covariates (linear regression). Differential Transcript Usage (DTU) analysis was performed using the R library IsoformSwitchAnalyzeR  (Vitting-Seerup & Sandelin, 2019)   (Anders et al. , 2012) ). For that reason, we chose isoformSwitchAnalysisPart1 function using DRIMSeq  (Soneson et al. , 2016)  to test differential transcript usage. Using part1 fasta files, all external analysis was performed and used as input to isoformSwitchAnalysisPart2 function. We used CPC2, Pfam, SignalIP and Netsurfp2 as indicated in the pipeline. Next, we performed a confirmation stage using stageR (Van den  Berge et al. , 2017)  to generate isoforms overall false discovery rate (OFDR). We selected all isoforms with OFDR < 0. 01 and dIF > 0. 05. RIN and PMI metadata were used as covariates (linear regression).
Single-cell RNAseq Using R library seurat, we created a seurat object (CreateSeuratObject), normalized data (NormalizeData), found variable genes (FindVariableFeatures) and rescaled data using a linear model (ScaleData, use. umi  = F). After that, we generated 50 PC's (RunPCA) but only used 35 of them based on the PC's visualization distribution (ElbowPlot). Since Allen data was already annotated, we only used tSNE (RunTSNE) to facilitate visualization. A group classified as "None" by Allen metadata were removed from our analysis. This strategy generated 7 main different cell types: Astrocytes, Endothelial cells, Glutamatergic Neurons, GABAergic Neurons, Microglia, Oligodendrocytes and oligodendrocyte precursor cells (OPCs). To assign genes to specific cell-types, we used the AverageExpression function.
RNAseq datasets obtained from different brain regions were used (Mayo: Allen et al, 2016; MSBB: Wang et al, 2018; ROSMAP: De Jager et al, 2018). Datasets were downloaded from AMP-AD Knowledge Portal (https://www. synapse. org) following all terms and conditions to use the data. The brain area analyzed and the number of individuals per condition were the following: Mayo -Temporal cortex, which neuroanatomically subdivides into the inferior, middle, and superior temporal gyri (STG), and cytoarchitectonically can be subdivided into Brodmann areas (BM, instead of BA) 20/21/22/41/42 (Strotzer, 2009), N=160 subjects (82 AD and 78 controls); MSBB -BM22, which is part of the Wernicke's area in the STG, N=159 subjects (98 AD and 61 controls); MSBB BM36, corresponding to the lateral perirhinal cortex, N=154 subjects (88 AD and 64 controls); MSBB BM10, corresponding to the anterior prefrontal cortex, N=176 subjects (105 AD and 71 controls); MSBB BM44, corresponding to the inferior frontal gyrus, N=153 subjects (90 AD and 63 controls); and ROSMAP -Dorsolateral prefrontal cortex (DLPFC), containing BM46 and part of BM9, N=423 subjects (222 AD and 201 controls).
Acknowledgments We thank  Charles Duyckaerts  (Neuroceb,  GIE Neuro-CEB   BB-0033-00011 ) for providing the brain samples. We would also like to thank  Dr. Tarciso Velho ,  Pierre Dourlen  and  Fabien Delahaye  for reading and suggestions on the manuscript. The results published here are in whole or in part based on data obtained from the  AMP-AD Knowledge Portal  (https://adknowledgeportal. synapse. org/). Study data were provided by the following sources: The  Mayo Clinic Alzheimer's Disease Genetic Studies , led by  Dr. Nilufer Taner  and  Dr. Steven G. Younkin ,  Mayo Clinic, Jacksonville , FL using samples from the  Mayo Clinic Study of Aging , the  Mayo Clinic Alzheimer's 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.
This work was also funded by the  Lille Métropole Communauté Urbaine  and the  French government's LABEX DISTALZ program  (Development of innovative strategies for a transdisciplinary approach to Alzheimer's disease) to JCL. DMC and LICC are supported by  Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES)  scholarships.
Data availability MSBB RNAseq (De Jager et al, 2018) -https://doi. org/10. 7303/syn3388564 ROSMAP RNAseq (Wang et al, 2018) -https://doi. org/10. 7303/syn3157743 MAYO RNAseq (Allen et al, 2018) -https://doi. org/10. 7303/syn5550404 Code availability Codes used for analyses are available, differentially expressed genes and differential transcript usage in individual datasets can be assessed using https://diegomscoelho. github. io/AD-IsoformSwitch/.
| Datasets and their assignment to groups We compiled 15 datasets describing invertebrate assemblages, flow permanence and human impacts in 406 rivers in eight European countries (Figure  2 ; see Appendix S4, Table  S1 ). We analysed all datasets in one all-region dataset, which identified interactions between aridity, impacts and river drying, and we therefore also analysed three regional datasets (based on biogeographical classifications;  EEA, 2002 EEA, , 2009 ; Figure  S1 ). We used nonmetric multidimensional scaling ordinations to visualize assemblage composition (Figure  S2 ) and to define coherent regional groups, removing one dataset characterized by exceptional compositional variability. Three datasets  (GB_NE, GB_SE and GB_SW)   S2c ). Three datasets were not assigned to regional groups due to insufficient replication and/or impact gradient length (Figure  S1 ; Table  S1 ).
We also analysed four individual datasets (CY, ES_NE1, ES_NE2 and ES_S), which represented all datasets with sufficient impact and drying gradients to warrant individual analysis  (Feld et al. , 2016; Figure   S1 ). All samples were collected during wet phases using standard quantitative or semi-quantitative methods designed to capture the taxa present in all habitat types (Table  S1 ). We excluded samples collected during non-flowing wet phases because flow cessation alters community composition  (Chadd et al. , 2017; Appendix S1. 1). We retained flowing-phase samples from all months to characterize responses to variability in environmental conditions across time, and explored seasonal differences in supporting analyses (Appendix S1. 2). We retained samples from perennial sites to encompass a full drying frequency gradient (0-1), and validated this approach in supporting analyses (Appendix S1.
3). We harmonized biological datasets to family level, which all datasets achieved for most taxa. We excluded Oligochaeta, which were typically identified to this coarser resolution, and excluded meiofauna and semi-aquatic taxa, which were recorded inconsistently. Family-level data can effectively characterize taxonomic and functional responses to human impacts and drying  (Datry et al. , 2014; Gayraud et al. , 2003) , and we validated this approach by comparison with genus-level analyses (Appendix S1. 4). We converted abundance data to presence-absence, to facilitate comparison of samples collected using different field methods  (Gayraud et al. , 2003 ; Table  S1 ), and validated this approach in supporting analyses (Appendix S1. 5).
| Functional trait assignment To identify traits promoting resistance and/or resilience to drying, we considered biological response traits in  Tachet et al. (2010) , plus female wing size and type, which describe the dispersal potential of aerial adult insects which have aquatic juveniles  (Sarremejane et al. , 2020). Informed by literature and expert opinion, we selected 22 trait categories across nine traits as responsive to drying F I G U R E 2 Study areas characterized by the 14 datasets used in all-region analyses, three regional analyses (bold text) and four individual analyses (dotted lines); freeform lines enclose larger study areas; CZ and HU datasets are nationwide. Datasets are labelled using ISO two-letter country codes, with abbreviated cardinal points describing the location relative to other national datasets (Table  S1 ).
1-Aridity index is the inverted ratio of precipitation to potential evapotranspiration; values increase with aridity (Tables  S4-S5 ). We weighted each category from 1 to 4, with higher scores indicating greater resistance and/or resilience to drying (Table  S4 ). For each genus, traits were coded using a fuzzy approach  (Chevenet et al. , 1994) , in which affinities between 0 (no affinity) and 4 (strong affinity) were assigned to each trait category. Affinities were then converted to percentage affinities within each trait and averaged across all genera with each family. We calculated family-specific resistance/resilience (RR) scores as the sum of each trait affinity multiplied by the corresponding trait category weight (Table  S4 ). Family scores were then ranked to identify the high RR group for each dataset, comprising its top 33% of families (Table  S6 ).
, 2002)  in the West Mediterranean region and STAR-ICMi  (Buffagni et al. , 2006)  in the East Mediterranean region. We calculated two functional metrics, redundancy (FuncRed) and richness (FuncRich), as described in Appendix S2. In brief, FuncRed was calculated as the difference between taxonomic diversity and functional diversity and represents the extent to which an assemblage is 'saturated' by taxa with comparable traits  (de Bello et al. , 2007). FuncRich was calculated as the multidimensional trait space representing each assemblage  (Villéger et al. , 2008). We also calculated functional metrics and taxonomic richness at genus level for selected datasets (Atlantic, ES_E, ES_S), and calculated FuncRed and FuncRich at a mixed subfamily level for the all-region dataset (Appendix S1. 4).