Shared brain transcriptomic signature in TDP-43 type A FTLD patients with or without GRN mutations.

Frontotemporal lobar degeneration with TDP-43 inclusions (FTLD-TDP) is a complex heterogeneous neurodegenerative disorder for which mechanisms are poorly understood. To explore transcriptional changes underlying FTLD-TDP, we performed RNA-sequencing on 66 genetically unexplained FTLD-TDP patients, 24 FTLD-TDP patients with GRN mutations and 24 control participants. Using principal component analysis, hierarchical clustering, differential expression and coexpression network analyses, we showed that GRN mutation carriers and FTLD-TDP-A patients without a known mutation shared a common transcriptional signature that is independent of GRN loss-of-function. After combining both groups, differential expression as compared to the control group and coexpression analyses revealed alteration of processes related to immune response, synaptic transmission, RNA metabolism, angiogenesis and vesicle-mediated transport. Deconvolution of the data highlighted strong cellular alterations that were similar in FTLD-TDP-A and GRN mutation carriers with NSF as a potentially important player in both groups. We propose several potentially druggable pathways such as the GABAergic, GDNF and sphingolipid pathways. Our findings underline new disease mechanisms and strongly suggest that affected pathways in GRN mutation carriers extend beyond GRN and contribute to genetically unexplained forms of FTLD-TDP-A.

Introduction

Frontotemporal lobar degeneration (FTLD), the second most common early onset neurodegenerative disorder, is highly heterogeneous. A definite diagnosis of FTLD can only be obtained after brain autopsy. The main pathology observed in FTLD patients is FTLD-TDP, named after the TDP-43 protein aggregates found in the brain. According to the type and localization of TDP-43 aggregates in the cortical layers, FTLD-TDP is classified into five pathological subtypes A to E, with A to C being the most common.[1,2] FTLD-TDP subtype A (FTLD-TDP-A) is characterized by abundant short dystrophic neurites, oval neuronal cytoplasmic inclusions in layer II of the neocortex and scattered lentiform intranuclear inclusions; FTLD-TDP subtype B (FTLD-TDP-B) patients on the other hand have few short dendritic neurites and moderate amount of neuronal cytoplasmic inclusions in deep cortical laminae; and in FTLD-TDP subtype C (FTLD-TDP-C) numerous long dystrophic neurites in upper cortical layers are reported, with low amount of neuronal cytoplasmic inclusion.

Several genes and genetic risk factors of FTLD have been described.[3-10] Importantly, there is a correlation between monogenic causes and FTLD-TDP subtypes. For instance, patients with a repeat expansion in the C9orf72 gene present mostly with FTLD-TDP-B at autopsy, while even more strikingly, patients with mutations in the gene encoding progranulin (GRN) invariably present with FTLD-TDP-A pathology.[11] Most FTLD-TDP-A patients, however, do not carry GRN mutations and are genetically unexplained. They represent ∼32% of the FTLD-TDP population.[8] We reported that homozygous carriers of the rare allele of rs5848, located in a microRNA binding site within the GRN 3′UTR,[12] are enriched in a genetically unexplained FTLD-TDP-A population as compared to controls with a 5-fold increased odds ratio.[8] The rs5848 risk allele is associated with reduced GRN expression in cerebellum, plasma and CSF.[12-14] Such correlations suggest a common aetiology between GRN mutation carriers and genetically unexplained FTLD-TDP-A patients. On the contrary, no genetic factors or specific pathways, have been associated with FTLD-TDP type C patients (FTLD-TDP-C).

Apart from genetic studies, important clues about the disease pathophysiology can be obtained from transcriptomic studies. Recently, a bulk RNA-sequencing study showed that Golgi vesicular transport, GABAergic signalling and cell death are affected in patients carrying C9orf72 repeat expansions as compared to non-expansion carriers and controls.[15] Also, Chen-Plotkin et al.[16] reported impairment of synaptic transmission and calcium signalling in GRN mutation carriers using a microarray approach.

To further dissect the transcriptional events underlying the different FTLD-TDP subtypes, with particular emphasis on the relationship between GRN-positive and GRN-negative FTLD-TDP-A, we performed bulk RNA-sequencing on the frontal cortex and cerebellum (a relatively less affected brain region) in a large cohort of FTLD-TDP patients and controls. We used multiple analytical methods including differential expression and coexpression analyses to identify genes and pathways associated with each FTLD-TDP subtype. We also aimed to identify targetable pathways relevant for FTLD-TDP patients.

Materials and methods

Participants

Participants with the frontal cortex and cerebellum available were selected from the Mayo Clinic Florida Brain Bank (n = 114; Table 1). Frontal-cortex tissue was collected from the middle frontal gyrus at the level of the nucleus accumbens. We included 24 FTLD-TDP-A, 20 FTLD-TDP-B and 22 FTLD-TDP-C, 24 GRN mutation carriers and 24 neuropathologically normal participants without any neurological disease (referred to as controls). FTLD-TDP-A patients had a median age at death of 83 [interquartile range (IQR): 79.0–86.8] and 12% were female. FTLD-TDP-B patients had a median age at death of 66.5 (IQR: 58.5–72.0) and 10% were female. FTLD-TDP-C patients had a median age at death of 74.5 (IQR: 67.8–78.0) and 8% were female. GRN mutation carriers had a median age at death of 66.5 (IQR: 64.0–75.8) and 14% were female. Controls had a median age at death of 86.5 (IQR: 79.3–89.8) and 16% were female. Patients with FTLD-TDP-A were all negative for GRN mutations as determined by Sanger sequencing. All patients were negative for C9orf72 repeat expansion using our previously reported protocol.[6] Cerebellum tissue was collected for a subset of these samples including 22 FTLD-TDP-A, 18 FTLD-TDP-B, 19 FTLD-TDP-C, 20 GRN mutation carriers and 21 controls. The tissues sampled did not include primary motor cortex. All GRN mutation carriers had type A histology at autopsy. Using principal component analysis on available genomic data, we determined that all samples were non-Hispanic white, except for two controls from African–American ancestry (Supplementary Fig. 1). The pathological diagnosis of FTLD-TDP was considered sufficient for inclusion, irrespective of the clinical diagnosis of the patient. However, 50% of FTLD-TDP-A patients were referred to the brain bank as clinical Alzheimer's disease, 21% as frontotemporal dementia (FTD) and 29% with other diagnoses that are part of the FTLD spectrum. GRN mutation carriers had a variety of clinical diagnoses with 20% referred to the brain bank as Alzheimer's disease, 20% with nfvPPA, 20% with FTD and 40% with other diagnoses that are part of the FTLD spectrum. Twenty-five per cent of patients with FTLD-TDP-B were diagnosed with behavioural variant FTD (bvFTD) with motor neuron disease, 20% with motor neuron disease alone and 60% with other diagnoses that are part of the FTLD spectrum. Thirty-six per cent of FTLD-TDP-C patients were diagnosed with bvFTD, 22% with svPPA and 42% with other diagnoses that are part of the FTLD spectrum. Genotypes of the rs5848 polymorphism were obtained by either direct Sanger sequencing of the GRN gene, or by a standard protocol TaqMan reaction with assay ID C_7452046_20 run on a QuantStudio 7 PCR System.

Table 1

Cohort characteristics

	Median age at death (IQR)	n female (%)	RIN FCX (IQR)	n total
Controls	86.5 (79.3–89.8)	16 (66.7)	9.1 (8.7–9.8)	24
FTLD-TDP-A	83 (79–86.8)	12 (50.0)	9.5 (9.0–9.7)	25
FTLD-TDP-B	66.5 (58.5–72)	10 (50.0)	9.8 (9.0–9.9)	20
FTLD-TDP-C	74.5 (67.8–78)	8 (36.4)	9.3 (8.7–9.7)	22
GRN mutation carriers	66.5 (64–75.8)	14 (56.0)	8.9 (8.2–9.5)	24

Information is shown for FTLD-TDP patients without mutation in the known genes, control participants without neurological diseases (Controls) and patients carrying a pathogenic mutation in the GRN gene (GRN mutation carriers). Age at death, RNA integrity number (RIN), number (n) of female participants and pathological diagnosis (FTLD-TDP-A, FTLD-TDP-B, FTLD-TDP-C) are specified. Data presented are median age at death [interquartile range (IQR)] or n and percentages (%).

RNA-sequencing

RNA from frontal-cortex tissue of 90 FTLD-TDP patients and 24 controls as well as from cerebellum tissue of 79 FTLD-TDP patients and 21 controls were extracted using the RNeasy Plus mini kit (Qiagen, Venlo, Netherlands). RNA quality and quantity were assessed using an Agilent 2100 Bioanalyzer and the RNA Nano Chip (Agilent Technologies). Only samples with an RNA integrity number (RIN) above seven were included in the study, with a median of 9.3 for the frontal cortex (IQR: 8.70–9.73) and 9.15 for cerebellum (IQR: 8.43–9.70). Library preparation was performed using Illumina TruSeq mRNA v.2 prep and sequenced at 10 samples/lane as paired-end 101 base pair reads on the HiSeq4000 (Illumina). Raw RNA-sequencing reads were aligned to the human reference genome (GRCh38) using the spliced transcripts alignment to a reference (STAR, v.2.5.2b).[17] Library quality assessment was performed using the RSeQC (v.3.0.0) package.[18] Gene-level expression was quantified using the featureCounts command in the Subread package (v.1.5.1).[19]

Differential expression analyses

Statistical analyses and plots were generated using open source R software packages available from CRAN and Bioconductor. Counts per million were obtained using the edgeR package using default parameters.[20] Samples were removed if their sex estimated from the RNA-sequencing data did not match the phenotypic information. Transcripts with low expression were removed using the filterByExpr function from edgeR package. Principal component and hierarchical cluster analysis using the Pearson correlation distance with average linkage were performed using the stats (v.3.6.2) and gplots packages. A negative binomial model was used to determine differentially expressed genes between groups. All analyses were adjusted for sex, RIN, batch and age. A Benjamini–Hochberg false discovery rate (FDR) correction was used for multiple testing. Genes with an FDR below 5% were considered statistically significant. Gene ontology (GO) analysis was performed using the anRichment package. GO terms with an FDR < 0.05 were considered significant. The up- and downregulated genes were used independently as input for STRING protein–protein interaction (PPI) network analysis. Only physical networks were considered based on text-mining, experiments and databases. Interacting proteins with a high confidence of 0.7 were used as input for Metascape PPI network analyses using default parameters.[21] In addition to using the STRING database, Metascape also queries the Biogrid database. The resulting networks and subnetworks were obtained with the subset of proteins that form physical interactions with at least one other list member. Using the means of Molecular Complex Detection algorithm, we identified connected subnetwork components. Each subnetwork was further subjected to GO enrichment analyses and drug target finding using the Opentarget database.[22] Venn diagrams were generated with the ggvenn package, volcano plots were created using the EnhancedVolcano package.

Coexpression analyses

We used the R package weighted gene coexpression network analysis (WGCNA) to identify sets of highly correlated transcripts (modules),[23] using as input the residual expression values adjusted for the aforementioned covariates. Residuals were obtained after conditional quantile normalization using the cqn package and the use of a multivariable linear regression model. Signed hybrid networks were created using the bi-weight mid-correlation method. To achieve a scale-free topology, we selected a power appropriate for each comparison (between 9 and 16). A dynamic tree cutting method was used with a minimum module size of 30 and a merged height varying from 0.2 to 0.3, depending on the comparison. Modules generated using these settings were represented by their first principal component (module eigengene) and a unique colour. For every gene, correlations between expression levels and each module's eigengene value (module membership) were calculated.

Modules that correlated with disease status were visualized using the top 20 most connected genes (hubs) as input for Cytoscape. In these networks, the connectivity of each gene was represented by the size of its node, the strength of the correlation is reflected by the thickness of its edges and their differential expression status is represented by the colour of the outer circle. The entire modules were then subject to enrichment analyses, as described before. Additionally, GO terms were collapsed using default parameters of the REViGO algorithm and the CirGO python library.[24,25] The top 500 hub genes from the green, purple and turquoise modules were used as input in the Metascape website for PPI analysis, as previously described. Networks were visualized with Cytoscape software,[26] followed by GO enrichment and drug target finding using the Opentarget database.

To construct consensus networks, all frontal-cortex data or only transcripts in common between frontal-cortex and cerebellum datasets were included, depending on the analysis. Similar parameters as the ones described previously were used for scale-free topology, module size and merged height. Correlation with disease status was performed as previously described. Fisher's exact test was used to compare the number of modules correlating with disease status in the FTLD-TDP-A and GRN consensus as compared to FTLD-TDP-C and GRN consensus. Conservation between modules across both tissues was assessed using a hypergeometric test for each of the pairwise overlaps.

Cell proportion estimation

Cell proportions were estimated using the R DSA package.[27] Marker genes were obtained using methods described elsewhere.[28] Briefly, we used the top 500 human specific markers from BRETIGEA R package for each cell type (neurons, oligodendrocytes, astrocytes and endothelial cells).[29] From this list, we randomly selected 100 genes for each cell type, and estimated cell proportions using the DSA algorithm. We repeated this 100 times and computed a Pearson correlation of estimated cell proportions between different runs. We then selected a representative subset of 100 marker genes per cell type as marker input for DSA analysis, which consists of computing the average marker gene expression in the dataset per individual. Differences in cell-type proportion per group of patients was assessed using a FDR corrected Kruskal–Wallis test.

We then applied the PSEA[30] method that allows model selection of the cell type(s) that should be included in controls or FTLD patients for each gene. Differential expression was subsequently estimated in specific cell types. For this, we used the functions em_quantvg and lmfitst and FDR correction. Finally, dysregulated genes for each cell type with a P-valueFDR < 0.05 were used as input into the IPA software for canonical pathway analysis using default settings (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuitypathway-analysis).[31]

Data availability

The datasets generated herein are available in the Synapse platform under accession number ‘syn25991921’. All other relevant data supporting key findings of this study are available in the article and Supplementary material.

Results

Differential expression of FTLD-TDP subgroups compared to control individuals

We performed RNA-sequencing on frontal-cortex tissue of FTLD-TDP patients (n = 66) without mutations in known FTLD-TDP genes including 24 FTLD-TDP-A, 20 FTLD-TDP-B and 22 FTLD-TDP-C, 24 patients with GRN mutations and 24 controls (Table 1). When performing quality control, two participants were excluded from this dataset based on sex discrepancies between reported and estimated sex, for a total of 88 patients and 24 controls.

Performing a principal component analysis, we found that the FTLD-TDP-A and GRN groups cluster together (Fig. 1A). Principal components 1 and 2 explained 34% and 9% of the variance, respectively (eigenvalue = 3.944 × 103, eigenvalue = 1.04 × 102, respectively). Confirming this finding, hierarchical clustering showed clustering of the FTLD-TDP-A and GRN mutation carriers (Fig. 1B).

Figure 1

Differential expression in FTLD-TDP groups versus controls. (A) Principal component (PC) 1 is represented on the x-axis and principal component 2 on the y-axis. Individuals with GRN mutation and pertaining to the FTLD-TDP-A group are coloured in pink and purple. Individuals pertaining to the FTLD-TDP-B and FTLD-TDP-C groups are coloured in light and dark green. Control individuals are coloured in blue. Ellipses delimit the different groups and are coloured accordingly. (B) Cluster analysis of all samples using the top 1000 most variable genes. Individuals with GRN mutation are coloured in pink, FTLD-TDP-A patients are coloured in purple. Individuals pertaining to the FTLD-TDP-B and FTLD-TDP-C groups are coloured in light and dark green, respectively. Control individuals are coloured in blue. Heat map rows show standardized expression levels of individual genes with red denoting high expression, yellow denoting medium expression and blue denoting low expression levels. (C) Volcano plots representing the differentially expressed genes in FTLD-TDP-A, GRN mutation carriers and FTLD-TDP-C groups versus control group. The fold change is presented in a log2 scale at the x-axis, while the FDR adjusted P-value is presented on the y-axis on a −log10 scale. (D) Venn diagram showing the overlap of up- and downregulated genes between groups (FTLD-TDP-A in purple, GRN mutation carriers in pink and FTLD-TDP-C in light green) as compared to controls. (E) Volcano plots representing the differentially expressed genes in GRN mutation carriers versus FTLD-TDP-C group.

When comparing FTLD-TDP-A to controls, and after adjustment for sex, RIN, age at death and plate, we detected 2704 differentially expressed genes (with P-valueFDR < 0.05 and at least two-fold change), from which 1139 were upregulated and 1565 downregulated (Fig. 1C and Supplementary Table 1). For the GRN group, 3380 genes were differentially expressed as compared to controls, from which 1221 were upregulated and 2159 downregulated (Fig. 1C and Supplementary Table 1). With these fold change and P-value thresholds, we detected only one differentially expressed gene in the FTLD-TDP-B group as compared to controls. We observed 791 differentially expressed genes in the FTLD-TDP-C group compared to controls, with 497 being upregulated and 294 downregulated (Fig. 1C and Supplementary Table 1).

Among the upregulated genes (found in the differential expression analyses of each patient group with controls), 31.0% were shared solely between GRN mutation carriers and FTLD-TDP-A compared to 6.7% shared between GRN carriers and FTLD-TDP-C (Fig. 1D): 15.8% of upregulated genes were shared between the three groups. Among the downregulated genes, 54.5% were shared solely between GRN mutation carriers and FTLD-TDP-A compared to 0.8% shared with FTLD-TDP-C (Fig. 1D): 9.3% of downregulated genes were shared between the three groups. AC104532.2 and AHNAK were the most upregulated and EFNA3 and SCN1B the most downregulated genes specifically in GRN mutation carriers and FTLD-TDP-A groups (Supplementary Table 1).

Since the presence of two risk alleles of the rs5848 variant, located in the 3′ UTR of GRN, is a known risk factor for FTLD-TDP-A, we next assessed whether the presence of homozygous rs5848 carriers was driving the similarities between GRN mutation carriers and FTLD-TDP-A. For this analysis, we removed all individuals carrying two copies of the rs5848 risk allele from all groups. This did not change substantially the principal component analysis and clustering results (Supplementary Fig. 2). In addition, when comparing groups of patients to controls, the strong similarity between GRN and FTLD-TDP-A was still obvious with 32.6% of the upregulated genes shared solely by GRN mutation carriers and FTLD-TDP-A, compared to 4.5% shared between GRN mutation carriers and FTLD-TDP-C (Supplementary Fig. 2). Among the downregulated genes, 40.8% were shared solely between GRN mutation carriers and FTLD-TDP-A, as compared to 1.7% shared with FTLD-TDP-C.

Even though we did not detect GRN differential expression in GRN mutation carriers as compared to controls in the frontal cortex, we detected an increased GRN expression in FTLD-TDP-A group (log2 foldchange = 0.48, P-valueFDR = 7.75 × 10−3). To assess differential expression in a brain region relatively unaffected by neuronal loss, we performed similar analyses in cerebellum tissue of patients and controls. We detected a reduced GRN expression in cerebellum in GRN mutation carriers although not significant (Supplementary Fig. 3, log2 foldchange = −0.22, P-valueFDR = 0.35).

Adding to our observations, when performing differential expression analysis on GRN mutation carriers versus FTLD-TDP-A patients, no significant genes were identified with P-valueFDR < 0.05 and log2 foldchange higher than 1 or smaller than −1. On the other hand, differential expression of GRN mutation carriers versus FTLD-TDP-C patients identified 1030 differentially expressed genes (Fig. 1E and Supplementary Table 1).

To further assess the similarities between FTLD-TDP-A and GRN mutation carriers, we performed consensus network analysis between GRN mutation carriers and FTLD-TDP-A or FTLD-TDP-C (Supplementary Fig. 4). Within the GRN and FTLD-TDP-A consensus analysis, we identified 21 modules, with eight correlating with disease status after Bonferroni correction (38.1%), while a consensus network between GRN mutation carriers and FTLD-TDP-C identified 23 modules, with only two correlating with disease status after Bonferroni correction (8.7%, P-valueFisher = 0.03). Altogether our results suggest that FTLD-TDP-A and GRN mutation carrier expression patterns are similar.

Differential expression analyses of combined FTLD-TDP-A with and without GRN mutations compared to control individuals

Given that our analyses suggested that FTLD-TDP-A and GRN mutation carriers share a common transcriptional signature, we combined the FTLD-TDP-A and GRN mutation carrier groups and compared them to controls. We identified 2902 differentially expressed genes (fold change < −2 or fold change > 2 and P-valueFDR < 0.05, Fig. 2A and Supplementary Table 2) from which 1235 were upregulated, including 766 genes that were part of the intersection between GRN mutation carriers and FTLD-TDP-A; and 1667 genes were downregulated including 1443 genes that were part of the intersection between the two groups. KANK2 was the most significant upregulated gene detected (log2 foldchange = 1.66, P-valueFDR  = 5.34 × 10−10). PCP4 was the most significant downregulated gene detected (log2 foldchange = −2.02, P-valueFDR = 1.76 × 10−11). Interestingly, we detected the two risk genes GFRA2 and DPP6 in the top 50 downregulated genes (log2 foldchange = −1.98, P-valueFDR = 2.02 × 10−10; log2 foldchange = −1.34, P-valueFDR = 2.82 × 10−10). We further investigated the most significant genes using enrichment analyses (Supplementary Table 3). Among the upregulated genes, we detected an enrichment for the terms ‘immune system process’ and ‘defence response’ (P-valueFDR = 2.74 × 10−45, enrichmentRatio = 2.30; P-valueFDR = 6.45 × 10−43, enrichmentRatio = 2.93, respectively). On the other hand, we detected among the downregulated genes an enrichment in ‘synaptic signalling’ and ‘nervous system process’ terms (P-valueFDR = 4.13 × 10−78, enrichmentRatio = 4.35; P-valueFDR = 3.80 × 10−38, enrichmentRatio = 2.96, respectively).

Figure 2

Differential expression analyses of FTLD-TDP-A with and without (A) Volcano plot representing the differentially expressed genes in FTLD-TDP-A and GRN mutation carriers combined, as compared to the control group. The fold change is presented in a log2 scale at the x-axis, while the FDR adjusted P-value is presented on the y-axis on a −log10 scale. (B) PPI subnetworks of genes negatively differentially expressed. (C) PPI subnetworks of genes positively differentially expressed.

We then performed PPI network analysis on the differentially expressed genes (up- and downregulated). Among the upregulated genes, we identified 10 subnetworks of interacting proteins (Supplementary Tables 4–6 and Supplementary Fig. 5). Of the 140 genes, 29 are known therapeutic targets for unrelated disorders. Among the downregulated genes, we identified 10 subnetworks of interacting proteins (Supplementary Tables 4–6 and Supplementary Fig. 6). Subnetwork 1 was enriched in ‘ion transmembrane transport’ terms (P-valueFDR = 6.93 × 10−5) and particularly in GABAergic genes (Fig. 2B). Out of the 203 genes, 87 are known therapeutic targets.

Coexpression analysis of combined FTLD-TDP-A with and without GRN mutations compared to control individuals

Next, we performed module-level analyses using WGCNA, adjusting for RIN, age at death, sex and plate. When comparing the combined group of FTLD-TDP-A and GRN carriers to controls, we identified 17 modules. Visualization of the module-trait relationships (Fig. 3A) revealed that the strongest relationships were dependent on the disease status with the identification of 11 modules of interest, after Bonferroni correction. As expected, no module correlated with potential confounders (RIN, age at death, sex or plate). GO enrichment analysis of these 11 modules (Supplementary Table 7) showed that they were involved in synaptic signalling (pink and turquoise, P-valueFDR = 4.90 × 10−3 and P-valueFDR = 8.18 × 10−58), Golgi vesicle transport (green-yellow, P-valueFDR = 2.67 × 10−2), RNA processing (lightcyan, P-valueFDR = 2.37 × 10−11), immune response (purple, P-valueFDR = 2.31 × 10−59), cardiovascular system development (green, P-valueFDR = 5.80 × 10−30), retrograde vesicle-mediated transport (salmon, P-valueFDR =1.87 × 10−3), RNA metabolic process (cyan, P-valueFDR = 1.30 × 10−7), SRP-dependent cotranslational protein targeting to membrane (blue, P-valueFDR = 2.23 × 10−31), and response to stimulus (yellow, P-valueFDR = 6.59 × 10−11). Among the genes in the green, purple and turquoise modules, a large number (43.85%, 45.55% and 44.56%, respectively, Supplementary Table 8) were differentially expressed and we therefore decided to focus the rest of the analyses on these modules. The eigengene of each module (Fig. 3B) showed that the green and purple modules were upregulated in patients compared to controls, while the turquoise module was downregulated in patients. More specifically, hierarchical clustering of GO terms revealed that the green module was enriched in terms associated with circulatory system development (41.8%) and immune response (21.4%, Fig. 4A). Although not a hub gene, GRN was upregulated (log2 foldchange = 0.33, P-valueFDR = 4.36 × 10−2) and belonged to this module. The purple module was enriched in GO terms associated exclusively with immune response (immune response 39%, positive regulation of immune system process 37%, Fig. 4B). Finally, the turquoise module was enriched in chemical synaptic transmission and regulation of ion transport GO terms (Fig. 4C). Interestingly, the hub gene for the turquoise module was NSF, which was identified in the PPI subnetwork 2 of downregulated genes (Fig. 2B). This PPI subnetwork was enriched in ‘vesicle-mediated transport’ terms (P-valueFDR = 6.23 × 10−3). To assess if the correlations observed in the frontal cortex were dependent on cellular composition, we performed consensus analysis on cerebellum tissue. Consensus network analysis between frontal-cortex data and cerebellum revealed that the green, purple and turquoise modules were conserved and correlated with disease status in cerebellum (Supplementary Figs 7–10).

Figure 3

Coexpression analyses of FTLD-TDP-A with and without (A) Module-trait relationships are presented for patients and controls. Modules with genes going up (red) or down (blue) together are shown, and a unique name has been attributed to each of these modules. Correlations and P-values are shown for variables of interest, including disease group (patients and controls), age at death, RNA integrity number (RIN), sex and plate. The strongest correlations (brightest colours) are observed for the disease group. (B) Heat maps are displayed for the green, purple and turquoise modules. High expression levels are shown in red and low levels in blue. Below every heat map, the first principal component of a given module (module eigengene) is displayed for each sample.

Figure 4

Top hub genes from the coexpression analyses of FTLD-TDP-A with and without (A) Top 20 hub genes for the green module are displayed. (B) Top 20 hub genes for the purple module are displayed. (C) Top 20 hub genes for the turquoise module are displayed. The size of the node is proportional to the connectivity of the gene. The strength of the correlation between two genes (nodes) is reflected by the thickness of the edge. Hierarchical GO analysis of biological process terms is presented on the right. Upregulated genes are circled in red and downregulated genes in blue. The main hub gene is coloured in yellow. Hierarchical GO analysis of biological process terms is presented on the right.

We then performed PPI network analysis on these modules of interest identified in the frontal-cortex WGCNA analysis. In the green module (Supplementary Fig. 11), we identified 13 subnetworks with the ITGA1 gene being part of subnetwork 4 and S1PR3 part of subnetwork 11 (Fig. 2C and Supplementary Tables 9 and 10). S1PR3 and ITGA1 were previously identified in subnetworks 1 and 4, respectively, from the differentially upregulated gene analysis. GO analysis revealed an enrichment of primary lysosome terms (P-valueFDR = 2.36 × 10−2) in subnetwork 8, which contains GRN (Supplementary Tables 9 and 10).

In the purple module (Supplementary Fig. 12 and Supplementary Tables 9 and 10), we identified nine subnetworks and more specifically, one enriched in major histocompatibility complex (MHC) protein complex (P-valueFDR = 4.59 × 10−2, subnetwork 1) and one enriched in neutrophil activation involved in immune response (P-valueFDR = 5.81 × 10−5, subnetwork 2). Interestingly, subnetwork 7 contains the C1q complex genes. Within the nine subnetworks, 14 out of 74 genes were reported drug targets (Supplementary Table 7).

In the turquoise module (Supplementary Fig. 13 and Supplementary Tables 9 and 10), we identified nine subnetworks and, similar to what we observed in the downregulated genes, we detected an enrichment in GO terms, such as regulation in calcium ion transmembrane transport (P-valueFDR = 1.48 × 10−3, subnetwork 2), and GABA receptor activity (P-valueFDR = 1.73 × 10−3, subnetwork 4). Subnetwork 4 contains the previously identified hub gene NSF. Within the nine subnetworks, 26 out of 125 genes were reported drug targets (Supplementary Tables 9 and 10).

Cell-type deconvolution of RNA-sequencing data

We further estimated the cell proportions in the frontal cortex of patients and controls. We observed that there was a strong variability between the different subtypes of FTLD patients, with a significant reduction of the relative proportion of neurons observed in GRN and FTLD-TDP-A patients, compared to controls (medianControls = 0.31, medianGRN = 0.04, P-valueGRN = 3.5 × 10−8; medianFTLD-TDP-A = 0.07, P-valueFTLD-TDP-A = 8.2 × 10−5; Fig. 5). We also observed a significant increase in the relative proportion of microglia in GRN compared to controls (medianControls = 0.03, medianGRN = 0.04, P-valueGRN = 2.7 × 10−2), and in the relative proportion of astrocytes in FTLD-TDP-A compared to controls (medianControls = 0.25, medianFTLD-TDP-A = 0.46, P-valueFTLD-TDP-A = 2.65 × 10−2). Finally, a significant increase in endothelial cell proportion was detected in GRN and FTLD-TDP-A compared to controls (medianControls = 0.06, medianGRN = 0.12, P-valueGRN = 5.10 × 10−4; medianFTLD-TDP-A = 0.12, P-valueFTLD-TDP-A = 4.34 × 10−3). These data indicate that the cellular composition in the frontal cortex of FTLD-TDP-A and GRN mutation carriers is similar, further highlighting the commonalities between the two groups.

Figure 5

Cell-type deconvolution analysis of FTLD-TDP-A with and without (A) Cell proportions estimations in the frontal cortex for five cell types (astrocytes, endothelial cells microglia, neurons and oligodendrocytes). (B) Residual expression in the frontal cortex and cerebellum of NSF per group.

Of the genes pertaining to the green module, 87% were predicted to be expressed by endothelial cells, of which 19% were specific for endothelial cells. Similarly, 97% of the genes pertaining to the purple module were predicted to be expressed by microglia, of which 18% were microglia-specific. Finally, 86% of the genes pertaining to the turquoise module were predicted to be expressed by neurons, of which 27% were neuron-specific.

We then performed differential expression in FTLD-TDP-A and GRN mutation carriers versus controls on the deconvoluted data. We detected 566 genes that were differentially expressed in neurons, 1144 genes in oligodendrocytes, 240 genes in microglia, 793 genes in astrocytes and 343 in endothelial cells (Supplementary Table 11).

Within neurons, PTGER3 and HECW1 were the most dysregulated genes (P-valueFDR = 1.59 × 10−9, beta = 0.08; P-valueFDR = 2.62 × 10−7, beta = −0.19, respectively) and the synaptogenesis signalling pathway was identified as the top perturbed canonical pathway (P-value = 2.05 × 10−14, Supplementary Table 11). Among other pathways, glutamate and GABA receptor signalling were also altered (P-value = 3.61 × 10−7 and P-value = 1.07 × 10−3, respectively). We detected two GABAergic genes in the top 10 differentially expressed genes (GABRB3 and GABRA5, P-valueFDR = 1.53 × 10−4, beta = 0.92; P-valueFDR = 1.53 × 10−4, beta = 0.33, respectively) confirming the implication of the GABAergic system in disease.

Within microglia, DDX46 and SLC22A9 were the most dysregulated genes (P-valueFDR = 4.15 × 10−4, beta = 0.86; P-valueFDR = 4.15 × 10−4, beta = −0.23, respectively). The GDNF family ligand-receptor interactions pathway was identified as the top perturbed canonical pathway perturbed (P-value = 1.58 × 10−3).

Within oligodendrocytes, GSTA1 and MAL were the most dysregulated genes (P-valueFDR = 1.73 × 10−5, beta = −0.03; P-valueFDR = 1.06 × 10−4, beta = −2.2, respectively). The most significant canonical pathway identified was the Ataxia-telangiectasia mutated (ATM) signalling pathway (P-value = 4.55 × 10−5).

Within astrocytes, WNT7B and TMCO4 were the most dysregulated genes (P-valueFDR = 1.65 × 10−4, beta = −0.17; P-valueFDR = 1.81 × 10−4, beta = 0.05, respectively). We detected an enrichment in mitochondrial dysfunction terms in differentially expressed genes (P-value = 5.24 × 10−6).

Within endothelial cells, HRCT1 and RNASEK-C17orf49 were the most dysregulated genes (P-valueFDR = 1.59 × 10−4, beta = −0.86; P-valueFDR = 4.71 × 10−4, beta = −2.15, respectively). The inhibition of matrix metalloproteases pathway was identified as the top perturbed canonical pathway (P-value = 1.72 × 10−4).

The NSF gene was identified in several of our analyses. It was significantly associated with disease status in the deconvoluted data in neurons (P-valueFDR = 2.15 × 10−2, beta = −0.79). Note that no significant difference in expression of NSF was observed in the deconvoluted data for FTLD-TDP-C and FTLD-TDP-B versus controls. We also highlighted NSF in the PPI analysis in the bulk differential expression analysis, as well as part of the top 20 hub genes from one of the three modules (green, turquoise or purple) and in their respective PPI analyses. The decreased expression of NSF in the bulk RNA-sequencing data appeared to be specific to the GRN mutation carriers and FTLD-TDP-A groups in the frontal cortex (Fig. 5B, log2 foldchange-FTLD-TDP-A = −1.66, P-valueFDR = 1.08 × 10−8; log2foldchange-GRN = −1.95, P-valueFDR = 2.78 × 10−7). We observed a trend in the same direction in cerebellum (log2 foldchange-FTLD-TDP-A = −0.31, P-valuenominal = 0.04; log2 foldchange-GRN = −0.38, P-valuenominal = 0.03).

Discussion

In this study, we used an unbiased whole transcriptomic approach to characterize expression patterns of FTLD-TDP patients in the affected frontal-cortex brain region. We investigated a total of 88 FTLD-TDP patients and controls. We found that GRN mutation carriers and FTLD-TDP-A patients have similar expression patterns using hierarchical clustering of the samples per gene expression and differential expression. We also reported that their frontal-cortex cellular composition is similar. Differential expression and coexpression network analyses identified an upregulation of inflammatory processes and downregulation of neuronal activity in patients compared to controls. Additionally, we highlighted vascular system, GABAergic and complement system alterations in FTLD-TDP-A patients. This study provides a comprehensive analysis of transcriptional changes in FTLD-TDP-A regardless of their underlying genetic aetiology and sheds light on potential new therapeutic targets and disease biomarkers.

GRN mutation carriers always present with FTLD-TDP-A at autopsy; however, it was not clear whether these patients shared common disease mechanisms with genetically unexplained FTLD-TDP-A patients. Previous genetic studies hinted towards a shared aetiology, reporting shared genetic associations such as TMEM106B and GFRA2 between both FTLD-TDP-A groups.[72] In this study, we showed that FTLD-TDP-A patients with and without GRN mutations shared a common transcriptional signature using unbiased approaches with principal component analysis and hierarchical clustering. We further showed more commonalities between GRN mutation carriers and FTLD-TDP-A than with FTLD-TDP-C patients regarding differentially expressed genes and coexpression networks. We previously reported that rs5848, located in the 3′ UTR of GRN, is a risk factor for FTLD-TDP-A, linking GRN to genetically unexplained FTLD-TDP-A[8]; however, we show that the similar transcriptomic profile was independent of rs5848. Furthermore, genetically unexplained FTLD-TDP-A patients have increased levels of GRN compared to controls, similar to findings in FTLD-TDP-C and FTLD-TDP-B patients, whereas GRN mutation carriers have levels comparable to controls. These data indicate increased levels of GRN in the frontal cortex as a result of activated microglia in all FTLD-TDP patients; however, since GRN mutation carriers have less GRN expression to begin with, their levels appear similar to controls.[32] These observations strongly suggest that the shared transcriptomic profile between FTLD-TDP-A and GRN mutation carriers is independent and likely downstream of loss of GRN. It is also noteworthy that we observed these highly similar transcriptomic signatures despite differences between FTLD-TDP-A and GRN mutation carriers with an average age at death 11 years later in the former as compared to the latter. Even though we observed overall similar coexpression patterns in cerebellum, a less affected brain region than the frontal cortex, the apparent similar transcriptomic signature is very likely to reflect similarities in cell-type composition. In fact, we observed a downregulation of PCP4 that encodes for PEP-19. PEP-19 binds to calmodulin, which activates several processes that are associated with neuronal death and play a role in the regulation of synaptic activity and maintenance of neuronal plasticity.[33] A downregulation of PCP4 may reflect the neuronal loss we observed in the deconvoluted data. Furthermore, deconvolution of our bulk RNA-sequencing data indicated a comparable cellular structure in the frontal cortex. In addition, the three highlighted modules were enriched in cell-type-specific genes suggesting that each cell type plays a role in disease. This observation is evocative of the cellular phase observed in Alzheimer's disease and argues strongly in favour of single-cell approaches to gain more granular insights into the pathophysiology.[34]

However, bulk RNA-sequencing is still a valuable resource to understand system level alterations. We highlighted the NSF gene that is significantly dysregulated in the frontal cortex of FTLD-TDP-A and GRN mutation carriers, even after deconvolution. NSF, which encodes for the N-Ethylmaleimide Sensitive Factor protein, is essential for vesicle-mediated transport through its interaction with the soluble NSF attachment protein receptor (SNARE) complex[35] and SNAPs.[36] The SNARE complex mediates synaptic vesicle fusion that is orchestrated by the assembly of three proteins: VAMP2, SNAP25 and STX1A and can be regulated by SNARE binding proteins such as complexin 1 (CPLX1).[37,38] In our data, we observed a downregulation of all these transcripts. Recently, an association between cognitive decline and dysregulation of the SNARE protein interactome, including CPLX1, was reported, with lower levels correlating with lower cognitive performance.[39,40] As CPLX1 is mainly expressed by GABAergic neurons, this indicates that impairment of the inhibitory system may contribute to cognitive decline.[41] To further elaborate, we also found a downregulation of GABAergic signalling pathways within the differential expression and the coexpression network analyses (turquoise module) and a dysregulation of GABA receptors in the deconvoluted neuronal population. In addition, the top dysregulated gene in the deconvoluted neurons, PTGER3, encoding for the Prostaglandin E Receptor 3, has been shown to act as a fever modulator through GABAergic signalling.[42] The GABA receptors expression patterns are highly complex and have been shown to present cell and regional specificity, suggesting that single-cell approaches would be appropriate to further dissect the GABergic system in FTLD.[43] A loss of GABAergic neurons has been reported in frontal and temporal cortices of FTLD patients, as well as reduced gamma oscillations and coherence in frontal lobes of bvFTD patients, reflecting lower GABAergic neuronal inhibition.[44,45] A recent study reported that GABA levels are reduced in the right inferior frontal gyrus of bvFTD.[46] While our observations are strongly associated with cellular composition, our data indicate that GABAergic synaptic vesicle fusion is impaired in FTLD-TDP-A patients. One could hypothesize that drugs targeting the GABAergic system might be beneficial.

ITGA1, that encodes the α1 subunit of integrin receptors, was identified in several of our analyses as upregulated in patients. Integrins are transmembrane heterodimeric receptors that are involved in cell adhesion and nerve regeneration.[47-49] The α1β1 complex is involved in neurite outgrowth through activation of the CDK5 signalling pathway.[50] Interestingly, the most upregulated gene identified in FTLD-TDP-A and GRN mutation carriers versus controls was KANK2. KANK2 has been shown to be the key protein involved in cell adhesion and migration through its interaction with integrins αβ.[51] In our deconvoluted data, ITGA1 was expressed by endothelial cells and part of the green module that was enriched in vasculature development terms. We also report an increased proportion of endothelial cells in FTLD-TDP-A and GRN mutation carrier groups as compared to controls, similarly to what is observed in Alzheimer's disease brains[52] that could explain the apparent upregulation of ITGA1. ITGA1 binds to type IV collagen[53] and COL4A2, which encodes one of the six subunits of type IV collagen, was also upregulated in our patient cohort. Type IV collagen is the major structural component of basement membranes and used as a marker for vasculature components in brain tissue. Intriguingly, Olofsson and Englund[54] recently reported, using this marker, a higher amount of a specific microvascular structure in the frontal cortex of FTLD patients. These raspberry-like structures are thought to be a sign of angiogenesis linked to hypoperfusion of the brain present in the bilateral frontal lobes of FTLD patients.[55-59] The green module was defined by a mixture between inflammatory signal and angiogenesis, suggesting a coregulation between both processes supporting the effect of microglial activation enhancing angiogenesis through pro-inflammatory TNF-alpha signalling.[60] Highlighting the importance of endothelial cells in FTLD, progranulin overexpression in endothelial cells impairs angiogenesis in vivo, increasing vessel growth and altering vessels’ integrity.[61] On the other hand, Grn knock-out mice have a profound disruption of the blood brain barrier associated with structural defects of the endothelial junctional complexes.[62] Among our deconvoluted data, we identified an enrichment in the inhibition of matrix metalloproteases pathway, with tissue inhibitors of metalloproteinases 4 and 2 (TIMP4 and TIMP2). Both play a critical role in extracellular matrix homeostasis, angiogenesis and inflammation.[63-65] Altogether, this suggests that angiogenesis is a fundamental part of FTLD-TDP type A pathology and further studies are necessary to elucidate its contribution.

Additionally, within the purple module, we highlight the C1q recognition component of C1, the classical complement pathway. C1q is a hexamer of trimers composed of A, B and C chains, encoded by the C1QA, C1QB and C1QC genes.[66] All three genes are upregulated in FTLD-TDP-A patients and GRN mutation carriers, probably reflecting increased numbers of microglia. Interestingly, C1QA has been shown to promote synaptic pruning by microglia in Grn−/− mice. In addition, knock-out of C1qa in Grn−/− mice protects against synapse loss, reduces behavioural deficits and improves survival, suggesting that complement activation and microglia-mediated synaptic pruning are drivers of neurodegeneration due to progranulin deficiency.[67] The C1q complex genes appear to be overexpressed in Grn−/− mice, and treating cultured Grn+/+ (wild-type) and Grn−/− neurons with C1q induces cytoplasmic TDP-43 granule formation.[68] Here we report that C1q upregulation is not only present in GRN mutation carriers but also in genetically unexplained FTLD-TDP-A patients.

Within the deconvoluted microglia, we highlighted the GDNF signalling pathway, which comprised the downregulated GFRA1 gene. GFRA1 is a RET coreceptor and part of the GDNF family ligands playing a neuroprotective role through modulation of microglial activation.[69-71] We previously reported a common variant at the GFRA2 locus as a genetic modifier for GRN mutation carriers,[72] expanding the number of GDNF family ligands involved in FTLD. In addition, GDNF implication appears to be common to several neurological disorders including Parkinson's disease where GDNF has been studied for its neuroprotective potential, even leading to clinical trials with GDNF brain delivery in patients.[73-75]

Within the deconvoluted oligodendrocyte data, we identified an alteration of the ATM signalling pathway. ATM, is a causal gene for Ataxia-telangiectasia (A-T) and encodes for a key protein involved in cell cycle checkpoint control during DNA damage repair.[76] The ATM kinase phosphorylates a number of downstream targets involved in DNA damage repair, cell cycle arrest and apoptosis.[77-80] Interestingly, A-T is characterized by cerebellar neurodegeneration and myelin abnormalities.[81] Recently, ATM deficiencies were shown to cause DNA damage accumulation in oligodendrocytes leading to a reduction of the mature oligodendrocyte lineage.[82] This is reminiscent of the TMEM106B−/− mouse model where a reduction of the number of differentiated oligodendrocytes was reported.[83] Along this line, TDP-43 plays a pivotal role in proper functioning of mature oligodendrocytes. Indeed, TDP-43 deletion in mature oligodendrocytes leads to progressive disruption of myelin and selective cell death of mature oligodendrocytes.[84] Interestingly, presence of TDP-43 inclusions in oligodendrocytes has been previously reported in FTLD-TDP, including GRN mutation carriers.[85,86] Altogether, our data add to the growing body of evidence that oligodendrocytes play a role in FTLD-TDP pathology.

Finally, the gene S1PR3, identified within several PPI analyses, encodes a G protein-coupled receptor for sphingosine 1-phosphate (S1P). Sphingolipids are highly abundant in the brain and play an essential role in neuronal plasticity.[87,88] We and others showed that GRN interacts with prosaposin,[89-91] a critical protein for sphingolipid metabolism.[92,93] In addition, GRN has been reported to bind β-hexosaminidase A and β-glucocerebrosidase, regulating their activities in sphingolipid metabolism and trafficking to lysosomes.[94,95] Importantly, S1P3 receptors are involved in neurite retraction and deceleration of axonal growth,[96] as well as cellular immunity function[97-99] and S1PR3 promotes inflammatory macrophage activation.[100] Knock-down of S1PR3 eliminates the anti-inflammatory effect of Fingolimod (FTY720), a neuroprotective drug used as immunosuppressant in multiple sclerosis.[101] Since S1PR3 is predominantly expressed by astrocytes, it was shown to be a crucial player in astrocytic inflammatory response on treatment with Fingolimod. Therefore, an increased expression of S1PR3 might increase neuroinflammation.

We also detected an alteration of RNA metabolism in the coexpression network analysis. This finding is in agreement with the dysregulation of spliceosome processes observed in neuronal cells on treatment with PGRN.[102] TDP-43 plays major roles in several aspects of RNA metabolism[103] and in particular, transcriptome analyses highlighted aberrant splicing in neurons with altered TDP-43 levels.[104,105] More recently alternative splicing was also reported in FTLD-TDP brains,[106-108] and additional studies are required to further dissect the specificity of splicing patterns observed in the different FTLD-TDP subtypes.

Our findings may have implications for other common neurodegenerative disorders. Indeed, TDP-43 pathology is present in late-onset Alzheimer's disease brains and is associated with memory loss and hippocampal volume.[109,110] Recently, a new pathological subtype of TDP-43 aggregates named TDP-43 type α has been described and consists of widespread distribution of dystrophic neurites and neuronal cytoplasmic inclusions in amygdala, hippocampus and frontotemporal cortex.[111] These features, reported in a subset of Alzheimer's disease patients, are similar to the TDP-43 immunoreactivity observed in FTLD-TDP type A.[111,112] Over half of these Alzheimer's disease patients were reported to have FTLD symptomatology. Further highlighting a potential similarity with FTLD-TDP-A subtype, the frequency of the homozygous TMEM106B protective allele (rs3173615) is reduced in non-FTLD brains with TDP-43 type α.[111] Our work now suggests that they may also share transcriptional profiles however, further studies are necessary to address this question.

Our study has some limitations. Even though this is a large cohort of FTLD-TDP patients studied by RNA-sequencing, the actual number of samples is still limited. Our findings of loss of neuronal activity and increased inflammation signalling might reflect end-stage molecular events. In contrast to our findings, Chen-Plotkin et al.[16] reported that GRN mutation carriers presented with a unique transcriptional signature. A possible explanation for this might lie in the higher number of brains studied here, and the use of more sensitive next-generation sequencing methods instead of microarrays. Additionally, we included only one affected brain region, whereas the aforementioned study used a mixture of two brain regions.

In sum, using an integrative approach from single gene differential expression to gene and protein network scale, we report for the first time a shared transcriptional signature in the frontal cortex of FTLD-TDP-A patients with and without GRN mutations. These transcriptional similarities are likely to be associated with comparable cell-type populations in both disease groups. We also describe an alteration in genes involved in vascularization through endothelial cells and importantly, we emphasize the dysregulation of the GABAergic signalling pathway, as well as the potential involvement of astrocytic inflammation. Overall, this study provides new avenues for research into FTLD-TDP type A and GRN and potentially new therapeutic targets.

Click here for additional data file.