Characterization of cell-cell communication in autistic brains with single-cell transcriptomes

Numbers and strengths of L-R interactions in autistic PFC and controls

Most current computational studies of cell-cell communications (or interactions) from single-cell transcriptomics data analyze ligand-receptor interactions among all cell types within one sample, but some software also provides methods for comparing interactions between two samples. We thus started with this common practice by comparing cell-cell interactions between ASD and control brains using a previously published snRNA-seq data [23] (Fig. 1a). The data contained samples from prefrontal (PFC) and anterior cingulate cortex (ACC), while cells/nuclei were classified into 17 types (see “Materials and Methods”). We focus our description on PFC below as its relation to social and cognition function and ASD is relatively better studied, but the same bioinformatics methods were applied to ACC data. We merged cells from all 13 ASD PFC samples into one group and cells from 10 control PFCs into the other and then studied the between-group difference (referred as “global approach”; Fig. 1a). We used the software CellChat, which analyzes the expression of ligands, receptors, and other proteins involved in cell-cell interactions to determine the interaction strength in an interacting protein complex for a pair of cell types [32]. For simplicity, we hereafter refer all pairwise interactions as ligand-receptor (L-R) interactions even though noncanonical ligands or receptors are also considered. In Additional file 1: Table S1a, we show the cell numbers in each of the 17 cell types and the numbers of L-R interactions in which a cell type acts as sender or receiver for both PFC and ACC. In total, CellChat identified 6433 and 5049 interactions in control and ASD PFC, respectively (Additional file 2: Table S2 a–b). The highest number of interactions in controls was 103, found in the L5/6-CC to L5/6-CC (cortico-cortical projection neurons) autocrine interactions, while the highest in ASD sample was 79, between L5/6-CC and L2/3 cell types. The total interactions in ACC are similar to those in PFC, with 5940 in control and 5843 in ASD (Additional file 2: Table S2 c–d). In both control and ASD ACC, L5/6-CC autocrine signaling had the most interactions, 126 in control and 109 in ASD. In both brain regions, L5/6-CC and L5/6 (layer 5/6 corticofugal projection neurons) are the cell types with the most interactions as either sender or receiver cells (maximal of 904 and 771 in L5/6-CC and L5/6, respectively), while microglia and Neu-NRGN-II have the least (maximal of 37 and 96 in microglia and Neu-NRGN-II, respectively). This difference is not directly related to cell numbers, because there is no significant correlation between cell numbers and L-R interaction numbers (Pearson’s correlation coefficient = 0.07) (Additional file 3: Fig. S1).

Although moderate (largest reduction is only 28, seen between L4 and L5/6 CC), we found a decreased trend (200 decreased vs 20 increased) in the numbers of cell-cell interactions in ASD vs controls for the 17 cell types in PFC (Fig. 1b; Additional file 3: Fig. S2a), involving both autocrine and paracrine signaling. This comparison of interaction numbers, however, has a caveat, because it could underestimate the contribution of strongly expressed ligands or receptors while amplify the small changes in lowly expressed genes; an interaction could be counted in one group but not in the other if its ligand or receptor is expressed in > 20% of cells (a threshold applied in CellChat to each cell type) in one group but < 20% in the other. The interaction strength, which captures the actual expression levels of all the genes encoding proteins in an interaction, is a better metric. With strengths of all L-R interactions considered (the maximal values are 0.37 in control and 0.4 in ASD), we found that cell-cell communications between many pairs of cell types became overall stronger in ASD (156 increased vs 133 decreased), with the increases most apparent for interactions involving excitatory and inhibitory neurons (Fig. 1c; Additional file 3: Fig. S2b). At the level of individual L-Rs, the interaction strengths also exhibited a trend of moderately increased in ASD (data not shown). The largest increase was found in an autocrine signaling for IN-SV2C interneurons, followed by signaling between IN-SV2C inhibitory neurons and L2/3 and L4 excitatory neurons, suggesting potential abnormal excitatory-inhibitory neuronal communications in ASD brains. Among the decreased interactions in ASD, paracrine signaling between oligodendrocyte and two inhibitory neurons (IN-SV2C and IN-PV) showed the most reduction, followed by OPC autocrine signaling. These findings are very interesting, not only in terms of major implications of inhibitory neurons but also with respect to the increasing appreciation of the role of glia in neuron development and functions, such as synapse modification [42, 43]. In addition, reduction of oligodendrocyte cells was previously reported and linked to neurodevelopmental disorders [44, 45].

The ACC shows a different pattern of alterations in cell-cell interactions, with similar numbers of cell-type pairs exhibiting increased or decreased interactions in ASD brains (Fig. 1d–e; also see below).

Major intercellular signaling patterns in PFC

Next, we grouped L-R interactions into signaling pathways using CellChat and compared them between ASD and controls. This is a key advantage of CellChat over other software, and it is important because many L-R interactions would invoke the same downstream signaling, such as FGF1/2/3-FGFR1/2/3 or NRNX1/2/3-NLGN signaling. Afterwards, we used CellChat to identify the major latent patterns in outgoing and incoming signaling in the control and ASD PFC, by clustering cell types sending or receiving similar signaling. For outgoing signaling (Fig. 2a and b), CellChat uncovered 5 and 4 patterns for control and ASD samples, respectively. In control (Fig. 2a), pattern 1 involved interneurons and Neu-mat cells, in which the main pathways were IGF, TAC, HGF, and OSTN. Pattern 2 involved solely excitatory layer neurons with NOTCH, WNT, EPHB, and NT as main contributors. Oligodendrocytes clustered by themselves in pattern 3, expressing SPP1, CD22, and MAG. Pattern 4 involved endothelia and microglia, with FN1, ESAM, OCLN, and PECAM1 as the coordinators. Finally, pattern 5 was comprised of only nonneuronal astrocytes and OPC cells, in which the expression of genes related to VEGF and ANGPTL was the most abundant. On the other hand, pattern 1 from ASD PFC was a mixture of nonneuronal and neuronal cells with astrocytes, interneurons, and Neu-mat (Fig. 2b). In this pattern, genes encoding IGF, CRH, BMP, and TAC signaling were the dominant ligands. In pattern 4, similarly, both nonneuronal and neuronal cells were involved, with a mixture of microglia and NRGN-expressing neurons. Here, CD39, CD45, and VISTA were the main senders. Like controls, excitatory neurons grouped together in pattern 2 with outgoing signals from NT, ANGPT, and CSF pathways. Pattern 3 clustered only nonneuronal cell types, endothelial and oligodendrocytes cells, with CD22, MAG, CLDN, and ESAM as the main signaling contributors.

Fig. 2figure 2

Cell-cell communication patterns in PFC. The networks show the CellChat inferred latent patterns connecting cell groups sharing similar signaling pathways. The thickness of the water flow represents the relative contribution of the cell group or signaling pathway to a latent pattern; outgoing patterns for secreting cells in control (a) and ASD (b); incoming patterns of receiving cells in control (c) and ASD (d). Note that CellChat only includes main contributors in these plots, and cell types (e.g., Neu-NRGN-I) making small contributions are thus omitted

From the receiving end, 7 patterns were found in control PFC while 6 in ASD (Fig. 2c and d). In control (Fig. 2c), pattern 1 involved excitatory layer neurons, while SEMA5, VIP, and WNT were the main signaling receivers. Pattern 2 grouped interneurons and Neu-mat cells with genes related to OSTN and PACAP pathways. FN1 and TENASCIN pathways were the main signal receivers in pattern 3, in which astrocytes participated. Endothelia, oligodendrocytes, and microglia clustered by themselves in patterns 4, 5, and 7. Finally, IN-SV2C were in pattern 6 with TAC, KIT, and CD226 as main receptors for incoming signals. For ASD (Fig. 2d), pattern 1 was formed by excitatory layer neurons with SEMA5, CCK, and SEMATOSTIN pathways as the main receptors. Interneurons and Neu-mat cells clustered again together in pattern 2, with TAC, KIT, and PACAP pathways. Genes involved in FN1 and TENASCIN pathways were the predominant receptors in pattern 3, enclosing astrocytes. Endothelial, both oligodendrocyte cell types and microglia were enclosed in patterns 4, 5, and 6, respectively, and showed same patterns as in control dataset.

Taken together, analysis of the intercellular CCC patterns indicates that all cell types contribute, although at various extents, in both control and ASD PFCs, and no cell type appears to serve as major signaling hubs in either condition. The CCC patterns seem similar in ASD and controls, but it appears a little less diverse (one fewer pattern) in ASD.

Pan-cell type signaling difference between ASD and control PFCs

Considering that one ligand may interact with receptors on multiple cell types, and vice versa, we next examined intercellular signaling for which multiple or all 17 cell types participated, in contrast to looking at cell-cell communications between any two cell types as described above. For this, we used CellChat to sum up all the L-R interaction strengths in the same signaling pathway across all cell types in which the L-R encoding genes were expressed. In total, CellChat has a collection of 229 families of signaling pathways. Among them, 68 unique ones were identified in our global approach. CellChat considers this pan-cell CCC network as “information flow.” Applying it to our data, we obtained “relative information flow” change between ASD and controls (Fig. 3a; Additional file 4: Table S3). The result indicates that 8 signaling pathways, noncanonical WNT (ncWNT), GRN, EGF, THBS, PTH, OCLN, SPP1, and CD226, were specifically identified in the control dataset, but no pathways unique to ASD were found. Additionally, 48 of the signaling pathways were significantly downregulated and 6 upregulated in ASD when the relative information flow was analyzed statistically (Fig. 3a; p < 0.05, Wilcoxon test). However, we should mention that the differences (i.e., effect sizes) for most signaling are relatively small (Additional file 4: Table S3). Close examination of the changes in relative contributions of individual cell types with respect to both outgoing and incoming signaling in ASD and control datasets illustrates this better (Fig. 3b and c). It also indicates that the two largest reduced signaling pathways (ncWNT and THBS) are due to significant decreases in incoming signaling in multiple cell types. This close-in analysis further suggests that excitatory and inhibitory neurons are the major sources of signaling senders, while both neurons and glial cells participate as receivers. The contribution of microglia to information flow, however, appears low in this analysis.

Fig. 3figure 3

Comparison of pan-cell type signaling networks in PFC. a Pan-cell type relative information flow showing signaling pathways identified in ASD PFC and controls. The pathways with greater information flow in ASD or controls were in cyan or red, respectively, with black indicating no significant differences. b and c Dot plots showing the difference in relative contribution of each cell type to outgoing (b) or incoming (c) signaling in ASD vs controls. Signaling on the left (a) with no difference was not included in (b) or (c)

In summary, our global analysis identified cell-cell communications and intercellular signaling that are potentially disrupted either between specific pairs of cell types or across many cell types in the ASD brains.

Sample-by-sample CCC analysis in PFC

The above global approach is standard in CellChat and is also the most common design in cell-cell communication analysis by other software (i.e., only making two sample comparisons). One concern is that it did not account for the potential differences among individual ASD or control samples. Furthermore, the computation would count interactions between cells in two different brains. Thus, we decided to apply a complementary sample-by-sample approach (Fig. 1). In essence, we used CellChat to obtain strengths of L-R interactions or score signaling pathways for each of the 13 ASD and 10 control PFC samples, as above, but then applied Wilcoxon test to determine if the sample-level strengths were significantly different between ASD and control groups. The numbers of cells for each ASD and control samples can be found in Additional file 1: Table S1b, while the pathway strengths for individual samples are in Additional file 5: Table S4a, including the mean values for ASD and controls, the standard deviation, and the differences as fold changes. We focused our analysis on pan-cell type signaling. In addition to the 68 identified by the above global approach, 24 others were detected. Of these 92 pathways, 10 were specific from controls and 3 unique to ASD. Pathways and genes belonging to each of them are listed in Additional file 5: Table S4b. After pathways found in less than 5 samples were excluded, 74 were retained for ASD vs control comparison. A principal component analysis (PCA) of the samples based on the scores of the 74 signaling showed that the first and second components explained 44.7% of the variation (Additional file 3: Fig. S3), with noticeable but not large separation between the ASD and control samples, indicating a small shifting but no major rewiring of intercellular signaling flow in the two groups, consistent with results above. To gain better insight into the underlying molecular interactions, we decided to analyze all the L-R interactions in the 74 signaling by Wilcoxon test. Although only 13 interactions remained significant at 10% FDR upon multiple testing correction, the results indicated that, at the nominal p < 0.05, 88 L-R interactions in these signaling were significantly different between ASD and controls (Fig. 4). These 88 L-R interactions were mapped to 32 signaling pathways (Additional file 5: Table S4c), 25 and 7 of which were down- and upregulated in the ASD PFCs, respectively, indicating that the other 42 signaling did not have significantly differential L-R interactions by our method. The data in Fig. 4 also help explain the apparent difference between an increase in pairwise cell-cell interaction strengths (Fig. 1c) and a reduction of global CCC information flow in ASD (Fig. 3a); there were more individual L-R interactions showing decreased strength (blue in Fig. 4), but they involved fewer cell types than those interactions exhibiting increased strength (red in Fig. 4).

Fig. 4figure 4

Heatmap for differential L-R interactions (y-axis) identified for individual pairs of cell types (x-axis) from the sample-by-sample approach

From the 32 differentially active pathways, we selected six with strong literature support for their potential involvements in ASD to illustrate our finding (see “Discussion” for details). They are neurexin (NRXN), fibroblast growth factor (FGF), contactin (CNTN), netrin-G ligand (NGL), neuregulin (NRG), and pleiotrophin (PTN) signaling. A chord diagram reflecting differences in interaction strengths between ASD and control samples is drawn for each pathway (Fig. 5), including the directions of changes and the cell types involved. For each of the pathways, we also included the genes if their expression levels were determined to be significantly different between ASD and controls cells by Velmeshev et al. [23]. Among the six, cell-cell interactions in NRXN and NRG were mostly increased in ASD while decreased for FGF and PTN, but similar numbers of increased and decreased interactions were observed in the other two pathways. Details of the strengths for each sample, cell pair, and pathway can be found in Additional file 6: Table S5.

Fig. 5figure 5

Chord diagrams plotting signaling strength differences between ASD and control PFC. The lines represent changes in L-R interaction strengths, with the statistically significantly different ones colored as intense red or blue, for increase or decrease in ASD, respectively. Light red or light blue for small changes not reaching statistical significance. Gray lines for no changes. Genes identified as differentially expressed in Velmeshev et al. [23] study were indicated in the corresponding cell type(s). The color bars in the inner circles indicates targeting cell types of the outgoing signaling while noncolor part for incoming signaling

To test how our result may be affected by the database of ligand/receptors, as it was reported to be a key factor [33], we replaced the CellChat L-R interactions with those from CellPhoneDB (v2.0) [31] and repeated the above ASD vs control comparison. This resulted in 43 signaling pathways, all of which were detected in the above analysis using CellChat database. Among them, 12 pathways showed a statistically significant difference between control and ASD (Wilcoxon test; p < 0.05).

In short, the results from our sample-by-sample approach strengthen our findings by the global method. Together, they highlight important intercellular signaling potentially disrupted in ASD PFC and demonstrate the advantage of CCC-focused analysis over conventional different gene expression analysis, as it can uncover altered intercellular signaling that would otherwise be missed in the latter.

Function enrichment analysis of the disrupted cell-cell signaling

To better understand the roles of the disrupted cell-cell signaling, we carried out GO enrichment analysis, using 165 genes in the 32 dysregulated pathways (Additional file 5: Table S4d). At FDR < 0.05, we found that the top “biological process”-related GO terms were axonogenesis, extracellular matrix organization, peptidyl-tyrosine phosphorylation, synapse organization, regulation of cell morphogenesis, among others (Fig. 6a) The enrichment analysis was reproduced by another software, ClueGO [35], which also reduced redundancy in the enriched terms. As shown in Fig. 6b, ClueGO was able to identify 8 clusters of GO terms, in which the main cluster, covering 27% of the terms, was related to axon development. This cluster is composed of neurogenesis, axon guidance, generation of neurons, neural crest migration or axon development terms, and pathways with important roles in neurodevelopmental disorders. The second most enriched cluster (24% of terms) was “enzyme linked receptor protein signaling pathway,” in which terms such as MAPK cascade, ERK1, and ERK2 cascade or protein tyrosine kinase activity were included. Regulation of neuron projection development covered 8% of the terms, comprising key pathways related to developmental disorders such as regulation of neurogenesis, positive regulation of axonogenesis, or regulation of nervous system development. Positive regulation of epithelial cell proliferation, biomineral tissue development, peptidyl-tyrosine phosphorylation, and vasculature development accounted for 9%, 5%, 4%, and 4% of terms, respectively. Taken together, the GO enrichment results are consistent with previous findings that have implicated these same cellular processes in ASD and other psychiatric conditions [10, 11, 13] but put a more clear connection to the important roles of cell-cell communication.

Fig. 6figure 6

Function enrichment analysis of genes in the dysregulated CCC signaling. a Dot plot showing the enriched GO terms. b Network connecting GO terms with sharing genes. Nodes are enriched GO terms, while the edges represent the extents of genes shared between two terms

Functional connection of the disrupted cell-cell signaling to ASD risk genes

In order to study how the disrupted cell-cell communication signaling is linked to ASD genetic risk, we applied the STRING database to construct a protein-protein interaction network. It identified 1076 edges connecting 153 of the 165 genes in the disrupted CCC signaling (ranging from 1 to 30 connections) to 239 ASD risk genes in the SFARI database (Fig. 7a). In addition, the network has an average local clustering coefficient of 0.38 and interaction enrichment values of 0.05, indicating overall strong connections among the genes (i.e., nodes).

Fig. 7figure 7

Connection between brain disorder risk genes and genes in the dysregulated CCC signaling. a Protein interaction network connecting dysregulated CCC signaling genes (blue) and ASD risk genes (white) in SFARI database. b Dot plot showing overlapping results of dysregulated signaling genes with lists of genes implicated in different brain disorders

Enrichment for gene sets related to brain disorders

To further address how disruption of the CCC signaling may be related to ASD, we intersected the 165 genes with lists of genes associated with ASD, schizophrenia (SCZ), intellectual disability (ID), bipolar disorder (BD), and attention-deficit hyperactivity disorder (ADHD) (Additional file 7: Table S6) and tested for significance of overrepresentation (Fig. 7b). Among the genes overlapping with ASD lists, CNTNAP2, LRRC4C, NLGN2, NLGN3, NRXN1, NRXN2, NRXN3, and PLXNA4 are considered as ASD risk genes with “high confidence” in SFARI database, due to multiple studies linked them to autism and other neurodevelopmental disorders [13, 46,47,48,49,50]. Also, significant enrichment was found with bipolar disorder (odds ratio (OR) = 2.7) and ADHD (OR = 2.4) gene lists, similar to previous findings [12, 51]. In addition to the ASD genes, SEMA5A or NRG1 are also very interesting. Semaphorin 5A gene has been linked to autism [52], and additionally, the failure of its expression has been linked to abnormal development of the axonal connections in the forebrain [53]. Decreased expression of neuregulin 1 has been linked to orchestration of a number of executive functions in ASD patients [54], as well as SCZ and bipolar disorders [55]. Regarding ADHD-enriched genes, Krumm et al. [56] related 3 semaphorin receptor protein (PLXNA3) to ASD. This receptor is known to be important for axon pathfinding in the developing nervous system [57]. Even though there were overlaps of genes in the ID and SCZ lists, no statistical significance was found for the overlap between the 165 CCC genes and genes related to these two disorders.

Relationship between disrupted intercellular signaling and differential gene expression programs

To address how disruption of the CCC signaling may lead to downstream gene expression difference in the ASD brains, we studied how they could be linked to the pathways enriched among DEGs in each of the cell types in the ASD brains. To do this, we first ranked genes by their expression fold changes between ASD and control cells in each of the 17 cell types and then performed GSEA (Additional file 8: Table S7). From the result, we selected and showed the enriched GO terms containing at least one gene in our disrupted CCC signaling (Fig. 8), resulting in 58 GO terms linked to 14 signaling pathways. Fourteen of the GO terms were more active in the controls, mostly involving nonneuronal cells. Oligodendrocytes showed 8 enriched pathways, including transmission of nerve impulse, regulation of glial cell differentiation, and cell chemotaxis. Similarly, AST-PP was enriched in integrin-mediated signaling pathways and regulation of smooth muscle cell differentiation. GO terms more active in ASD samples were mainly seen in neuronal cells, i.e., IN-SST, IN-SV2C, IN-VIP, L2/3, and Neu-NRGN-I. In total, 45 GO terms more active in ASD neurons were connected to disrupted CCC signaling, including synapse assembly, neuron cell adhesion, learning, regulation of neuron differentiation, postsynaptic specialization organization, receptor clustering, or synaptic signaling. Taken together, these results indicate that disrupted intercellular communications can be linked to downstream intracellular pathways and gene expression programs dysregulated in both neurons and glia cells in the ASD brains.

Fig. 8figure 8

Connection between cell-type ASD-enriched pathways and dysregulated CCC signaling. Dot plot showing enriched pathways from GSEA for individual cell types, with red and blue for higher activities in ASD and control PFC, respectively. The right column lists the corresponding dysregulated signaling from CCC analysis

Spatial co-expression of ligands and receptors

Our analysis assumed that all cells of the same type express the same ligands and receptors, by using snRNA-seq data that did not directly include spatial locations of cells, although the excitatory neurons did contain brain layer information. Spatial location, however, could be an important factor in cell-cell communication. To address this, we analyzed the spatial transcriptomics data from Maynard et al. [41], which applied 10× Genomics Visium platform to generate a map of gene expression in the six-layered dorsolateral prefrontal cortex of an adult human brain. In their study, the authors have already corroborated their spatial registration to the cell-type annotation in the snRNA-seq data used in our analysis. Here, we wanted to determine if the spatial expression levels of the ligand-receptor pairs are more correlated than expected by chance. The data indicate that most of the randomly paired L-R had correlation coefficients near “0,” but the interacting pairs identified by CellChat showed significantly larger correlations, using spatial data across all brains or the subset in defined layers (Additional file 3: Fig. S4a/b). The correlation shifted to both positive and negative perhaps because CellChat interaction includes soluble agonists, antagonist, and co-stimulatory, and co-inhibitory cofactors.

Cell-cell communications in anterior cingulate cortex

We present our findings for PFC above, but we have also carried out analyses for the anterior cingulate cortex samples using our global approach. Similar to PFC, most of the cell-cell interactions in ACC also occurred through neuronal cell types (e.g., L5/6-CC and L2/3 excitatory neurons), while microglia and Neu-NRGN-I contributed the least. In contrast to PFC, our global approach found a mixture of lost and gained pairwise interactions in the 17 ASD vs controls as mentioned above (Fig. 1d; Additional file 3: Figure S5a). L4, Nue-mat, Neu-NRGN-I, Neu-NRGN-II, and AST-PP cell types showed an increased number of interactions, acting as receiver cell types, while the remaining cell types showed decreased interactions in ASD. Similar to PFC, CCC strengths (Fig. 1e; Additional file 3: Fig. S5b) were increased in ASD ACC for most of the cell types. For pan-cell type signaling network, CellChat identified 70 unique signaling pathways in ACC, with 5 specific to controls (ANGPT, CD226, CD99, ENHO, and ncWNT) and 4 to ASD (EDN, FN1, PTH, and TENASCIN) (Additional file 9: Table S8). In Additional file 3: Fig. S6, we show the relative information flow for these pathways, with 37 significantly up- and 16 downregulated in the ASD. Only 13 pathways showed a decreased activity, while 5 showed increased activity in both ASD PFC and ACC. Further analysis compared the relative contributions from each cell type for outgoing and incoming signaling in ASD and control ACC, indicating global similarity but noticeable changes for some pathways, such as CD99, ncWNT, or ANGPTL (Additional file 3: Fig. S7a/b).

留言 (0)

沒有登入
gif