Differential Transcriptional Programs Reveal Modular Network Rearrangements Associated to Late-Onset Alzheimer’s Disease

1. Introduction

Alzheimer’s disease (AD) is a clinical-pathophysiological entity, currently conceptualized as a continuum [1]. It is marked by the presence of dementia symptoms alongside specific molecular biomarkers. AD is commonly categorized as early-onset (EOAD) or late-onset (LOAD) based on an age cut-off, typically 65 years [2]. LOAD is the most prevalent neurodegenerative disorder in older adults. It affects 5%–10% of this population, with prevalence rising to 50% in individuals over 85 years old [3]. Aging is the most important risk factor for the development of AD [4]. The rise in the aging population, along with enhanced longevity, underscores the imperative for the prevention and treatment of progressive neurodegenerative diseases, particularly LOAD [5].

LOAD has been linked to disruptions in multiple biological pathways, often occurring simultaneously during disease progression. The condition presents significant complexity in both its underlying causes and clinical symptoms. The diverse molecular phenotypes observed in complex diseases such as LOAD arise from the coaction between genetic predispositions and environmental factors, which together influence various pathways and processes. For example, genetic variants associated with LOAD can impact multiple biological systems, including immune dysfunction, altered lipid metabolism, and neuroinflammation [6,7]. Additionally, a spectrum of other prevalent and complex disorders has been linked to an increased risk of developing Alzheimer’s disease [8,9,10,11,12], including other psychiatric disorders [13,14,15]. Research in this area has significantly advanced our understanding of the biological and environmental foundation of Alzheimer’s disease, particularly its association with $\beta$–amyloid and tau-related mechanisms. However, the precise molecular events and biological pathways driving the disease remain elusive. This genetic heterogeneity poses substantial challenges for the development of targeted therapies and the identification of disease-associated molecular changes. Moreover, the limited understanding of how these pathological molecular events interconnect further impede the discovery of effective therapeutic targets. Given the intricate genetic interactions involved in LOAD, the condition presents a compelling case for investigation through network analysis.

Post-mortem diagnosis of AD is staged by the severity and distribution of these pathological hallmarks: extracellular $\beta$-amyloid deposits and intracellular neurofibrillary tangles (NFTs) in neurons [16]. Tangles initially emerge in the entorhinal cortex (Braak stages I–II), subsequently spreading to the hippocampus and thalamus (Braak stages III–IV), and ultimately reaching the neocortex (Braak stages V–VI). This progression generally correlates with the worsening of cognitive function, ranging from mild cognitive impairment to severe dementia [16,17,18]. The dorsolateral prefrontal cortex (DLPFC) is located in the middle frontal gyrus (MFG) of the frontal lobe of the neocortex. It’s an area in the prefrontal cortex, which is not an anatomical structure, but a functional one, located in the lateral part of the 9 and 46 Brodmann areas (BA) [19]. DLPFC does not operate in isolation; it interacts with various brain regions, including the thalamus [20], hippocampus [21], and associative areas like the posterior temporal and parietal cortices [22] and it is recognized for its role in working memory [23], speech [24], cognitive control [25], cognitive flexibility [26], planning [27], impulse inhibition and self-control, risk/benefit evaluation [28], problem solving [29], among other high order cognition functions. These functions are impaired during AD progression [30], and this region has been proposed as a therapeutic target to improve working memory in individuals with AD [31]. For these reasons, the DLPFC was selected as the focus of our analysis due to its role in cognitive processes most susceptible to impairment in AD and its potential as a target for therapeutic interventions.

1.1. Transcriptional Networks

Networks are basic tools in systems biology for expressing the components of biological systems as whole integrated systems [32]. Many natural systems are structured in the form of networks. For example, genes do not function in isolation, but rather as part of large and intricate information processing networks that involve multiple molecular entities [33,34] involving dynamic events, which are subject to regulation and integration at multiple levels, each with multiple steps. In this way, gene expression is a complex process that requires the coordinated regulation of multiple steps, from DNA to RNA, including transcription initiation, elongation, termination, and post-transcriptional modifications. To better understand this complexity, gene co-expression networks offer a powerful framework to investigate intricate biological systems, such as expression landscapes in disease [35].

Gene expression can be structured as a matrix thus characterizing gene expression data as multivariate. When converting multivariate data into a network, it becomes essential to apply metrics or some form of measurement capable of quantifying the distance or similarity between two models [36]. Using a correlation measure, we can examine how genes are co-expressed within an expression matrix, resulting in the construction of an adjacency matrix that is mathematically homomorphic to a co-expression network. Transcript co-expression networks often exhibit complex topological structures and a high degree of heterogeneity in their connection strengths and capacities. These networks are distinguished by the presence of “hubs”—nodes with a disproportionately large number of connections compared to the average degree and therefore have a robust architecture as they are resistant to random node failure [37]. However, the presence of hubs makes scale-free networks vulnerable to targeted insults where a hub is compromised [38]. In co-expression networks, hub genes frequently contribute to the stability and organization of the system. Furthermore, gene organization and expression are structured within modular networks [39], where certain genes within these communities of genes can serve as drivers of specific functionalities and organizational patterns. By analyzing relevant network metrics, the graph’s topological structure can be leveraged to gain insights into the macro, mesoscopic and local landscape of health and disease. These network relationships, when considered collectively, offer a snapshot of a phenotype’s state and provide a basis for comparing distinct phenotypes [40].

Transcriptome network analysis has been a tool in understanding large-scale dysregulations associated with complex diseases such as AD and dementia [41]. In addition, downstream network analysis has been useful in nominating potential therapeutic candidates and uncovering disease-specific patterns [42,43] and recognizing genes related to medical characteristics [44,45]. We aim to identify key metrics that characterize the disease-associated network structure and organization, ultimately highlighting the biology (genes, transcripts and biological pathways) involved in the pathological processes of late-onset Alzheimer’s disease (LOAD), specifically in the dorsolateral prefrontal cortex (DLPFC). Building on these observations, this study presents a statistical analysis of a co-expression network model of the DLPFC associated with Alzheimer’s disease. By leveraging RNA-seq data and probabilistic modeling, we develop both functional and semi-mechanistic models of the disease.

2. Methods

2.1. Data Acquisition and Classification

A complete workflow for this study is depicted in Figure 1. RNA-seq count matrices from the dorsolateral prefrontal cortex (DLPFC), specifically encompassing Brodmann area BM46 and part of BM9, were obtained from 793 subjects (473 with AD and 301 controls) through the Religious Orders Study and Memory and Aging Project (ROSMAP) database [46]. These datasets were downloaded from the AMP-AD Knowledge Portal (https://www.synapse.org) in compliance with all applicable terms and conditions.

RNA sequencing was performed following the methodology described in De Jager et al. [47]. Libraries were prepared using rRNA depletion, and sequencing was conducted on an Illumina NovaSeq 6000 platform. Metadata was used to classify patients into Control and AD groups according to the AD neuropathologic changes, following the NIA-AA criteria as outlined in the revised guidelines for the neuropathologic assessment of brain disorders common in the elderly, including Braak staging for assessing the severity of neurofibrillary tangles and the Consortium to Establish a Registry for Alzheimer Disease (CERAD) semi-quantitative measure for neuritic plaques (CERAD score) [48].

2.2. Quality Control

Quality control of the data was performed using the NOIseq package in R [49], 2015). After a diagnostic analysis, transcript composition bias, and length bias were corrected. Genes with low counts were filtered using a threshold of 1 CPM. Batch effects were corrected using the ASCA Removal of Systematic Noise on Seq data (ARSyNSeq) method, with the ’sequencingBatch’ column as the reference. ARSyNSeq filters the noise associated with identified or unidentified batch effects considering the experimental design and applying Principal Component Analysis (PCA) to the ANOVA parameters and residuals [50]. After this processing 28,263 features remained, including protein-coding genes, long non-coding RNAs, microRNAs, pseudogenes, and other RNA species. The entire dataset was then normalized by the trimmed mean of the M-values normalization method (TMM), which also adjusts by transcript length [51] (See Supplementary Figure S1). Finally, data was stratified by diagnosis, resulting in two expression matrices.

2.3. Inference Of Co-Expression Networks

Mutual Information-based gene coexpression networks were constructed separately for two groups: individuals diagnosed with Alzheimer’s Disease (AD) and age-matched controls without AD (control). Nodes in the network represent genes, and edges represent significant coexpression relationships between them based on Mutual Information (MI). MI is the maximum entropy/maximum likelihood estimate of statistical dependence between two random variables [52]. MI has attractive information-theoretic interpretations and can be used to measure non-linear associations in complex systems [53], such as transcriptional programs. This characteristic is advantageous since many relationships between transcripts do not follow a linear statistical dependence. As MI is well-defined for discrete or categorical variables, we first discretized the counts and posteriorly calculated the mutual information by gene pairs with the Infotheo R package [54]. To establish a MI heuristic cut, basic network metrics by percentile were obtained for percentiles 99.9999%, 99.999%, 99.99%, 99.9%, 99%, 98%, 90%, and 80% (See Supplementary Figure S2). Also, we ran a null model to certify that the network had characteristics not due to randomness ((See Supplementary Figure S3)

By analyzing various network metrics across different percentiles (e.g., number of vertices, edges, components, clustering coefficient) provides a clear picture of how the structure changes as we move from lower to higher MI interactions. The 99.99% percentile retains a manageable number of nodes and edges while maintaining the largest components and preserving important topological features while filtering out less meaningful interactions. By selecting interactions at the 99.99% percentile, we ensure that only the statistic dependence (i.e., those with the highest MI values) are included in the network. This helps to filter out noise and weaker associations that may not contribute meaningfully to the biological interpretation. This approach is equivalent to stratifying the data based on p-value thresholds [55]. The networks were then analyzed using the igraph package in R version 2.0.3 [56] and Cytoscape v3.10.2 [57].

2.4. Topological Analysis and Network Centralities Measure

A common approach for analyzing specific network properties is to compute and compare certain local or global network indices such as degree distribution and connectedness [58]. We began by identifying key quantities that characterize the structure and organization of the networks. First, we evaluated the degree distribution behavior in both models. The degree of a node represents the number of direct connections it has with other nodes. We further compared the cumulative degree distributions between the two models to assess topological differences. A Kolmogorov Smirnov (K-S) test was performed to test whether both distributions were significantly different. We performed a goodness-of-fit test to determine the best-fitting distribution for our data. Various distributions, including power-law, Poisson, normal, Pareto, exponential, and negative binomial, were evaluated based on statistical criteria (See Supplementary Tables 1 & 2)

To compare the overall connectivity differences between the models, the Jaccard similarity for the edges was calculated by comparing the sets of edges in each network. The computed index shows the ratio of the intersection of edges between the two networks to their union, as described in Equation 1.

$$J(E_{g_1},E_{g_2})={\displaystyle \frac{|E_{g_1}\cap E_{g_2}|}{|E_{g_1}\cup E_{g_2}|}}$$

(1)

Equation 1: The Jaccard Index equation was used to compare the edges between both networks.

Genes with high connectivity based on the degree and betweenness of nodes in the models were identified. In coexpression networks, Hub genes are genes with a significantly high number of direct connections to other genes, serving as hyper-connected nodes that contribute to network stability and functional coordination. We defined hub genes as those with the highest node degree, specifically the top 10% of the degree distribution. Moreover, betweenness centrality measures how often a node appears on the shortest path between pairs of other nodes, indicating its importance in controlling information flow across the network [38]. By occupying a central position in the network, they influence the flow of biological information, signals, and resources, ensuring efficient communication between different regions. Consequently, they can either facilitate interactions or act as bottlenecks within the system. Accordingly, high betweenness centrality genes were defined as those ranking in the top 10% of betweenness values. To explore the biological functions of the AD-specific hub and high betweenness genes, we performed a GO enrichment analysis using the clusterProfiler package in R. We focused on Biological Process (BP) ontology terms with a p-value cutoff of 0.05 and a q-value cutoff of 0.05. To gain a deeper understanding of connectivity changes related to hubs, we identified and enriched the first neighbors of the hubs specific to the AD network that were also exclusively first neighbors of hubs within the same network.

To visualize the connectivity of AD-specific hub genes, we extracted an induced subgraph comprising only the hub nodes from the AD network. This subgraph highlights the dense interactions between these key genes (see Figure 2c).n parallel, an induced subgraph comprising only the high betweenness nodes from the AD network was built (See Figure 2d).

Transitivity (an approximation to the the global clustering coefficient) [58] was computed for both control and problem networks. This is a measure used to detect the fraction of triplets that are closed in directed networks [59]. In other words, transitivity measures the probability that the adjacent vertices of a vertex are connected [60].

2.5. Inference of Modular Structure and Functional Analysis

Modularization aims to solve a complex problem by breaking it down into conceptually manageable parts. With the increasing scale and complexity of biochemical networks, a shift toward a top-down approach has occurred. In this framework, previously characterized reactions are integrated into modules based on their co-occurrence in wiring diagrams [61] associated with specific phenotypes. Such modules facilitate the study of biological processes by deconstructing complex biological networks into conceptually simple entities. To identify modular substructures in both models, networks were partitioned into gene modules using the Infomap algorithm for community detection [62]. Infomap is based on information theory and works by representing the movement of a random walker through the network to capture the modular structure in networks [63,64].

We used normalized mutual information (NMI) as a community detection comparison metric to assess differences in the modular structure of AD and control networks. NMI was computed using the compare() function in the igraph package. Then, to evaluate the similarity of genes and biological functions across the entire network based on module enrichment, we constructed two similarity matrices: one representing gene correspondence to modules and another capturing enriched modules according to their associated biological functions. Functional similarity was quantified using the Jaccard index of semantic similarity scores, allowing for an objective comparison of the biological processes represented within each module by their respective functions.

To determine the representation of each biological process across the network, we quantified how many modules (or communities) are associated with a specific biological process following Gene Ontology (GO) enrichment analysis. We identified the number of modules in which each biological process is represented, providing insight into the distribution of these processes across the network (See Figure 3).

3. Results

Both graphs showed a power-law-like behavior, signifying that a small number of nodes, or hubs, are highly connected, while the majority have fewer connections (See Figure 4). Both networks are essentially disassortative. The Scaling exponent ($\gamma$) in the log-log plot indicates that the AD network has a slightly higher slope, implying that the AD network tends to have more nodes with a higher number of connections than the control network. This is confirmed in the cumulative frecuency representation, where we see the same phenomenon but represented in a cumulative degree graph.

In the AD network, certain genes or nodes play a central role in maintaining connectivity across the system. In contrast, the control network shows a slightly more homogeneous and decentralized structure. This reflects a more balanced distribution of connectivity, where fewer nodes have a large number of connections, leading to a network that is less dependent on a small set of hubs for its overall cohesion. The K-S test showed that there is no statistically significant evidence to claim that the underlying distributions of the two samples are different (D = 0.037916, p-value = 0.9971). Table 1 shows the topological features of both networks.

Edges of the AD and control networks shared only 68.39% similarity. Also, two key centrality metrics, degree centrality and betweenness centrality, were analyzed to identify key genes in the topology of the networks. For both metrics, the 10% of genes with the highest degree and betweenness were taken. We found 112 Hub genes in the AD network and 108 in the control network. Both models share 82% of the Hub genes, but there are 11% of genes found only in the AD network. In parallel, 134 and 129 High betweenness centrality genes were found in the AD and control network, respectively. 52% of the high betweenness centrality genes are common to both conditions, while 25% of the genes were found only in the AD network (see Figure 2a,b). Only the high betweenness genes enriched pathways in the Biological Process (BP) ontology test. Pathways enriched included primarily the synaptic vesicle cycle, synaptic vesicle maturation, regulation of synaptic vesicle exocytosis and vesicle-mediated transport in synapse (see Figure 2e).

3.1. Mesoscopic Network Analysis

Graph partitions showed 68 modules of genes being coexpressed for the AD graph and 71 modules for the control graph (see Figure 5). By comparing the Q value for each network partition it is possible to evaluate whether there is a network in which the module detection algorithm, applied to each network, identifies a more modular structure [65]. The modularity of the AD graph is slightly higher (0.28) compared to the control graph (0.20), suggesting that the AD network has a clearer community division. We further analyzed the modular structure of the networks by comparing the community partitions derived from the Infomap algorithm, which yielded a NMI of 0.4479.

A similarity matrix based on the membership of genes to their respective modules was constructed to understand the rearrangement of genes into different modules. Of the modules in both networks, 10 modules were preserved with perfect similarity between the two networks. In addition, 11 modules showed a similarity of more than 90%, 12 modules had a similarity greater than 80%, 14 modules exceeded 70%, 18 modules showed a similarity greater than 60%, and 24 modules achieved a similarity of more than 50%. This points to a reorganization of genes in several modules between the two networks. Although some of the modules remain largely intact, in others a significant redistribution of genes between different communities is observed (See Figure 6a).

Other similarity matrix, now that compares the enriched biological functions of individual modules between the two networks to evaluate the similarity of biological functions across the networks (Figure 6b). This compares the enriched biological functions of individual modules between the two networks. Gene Ontology (GO) terms for Biological Processes (BP) were assigned to each module, allowing for a detailed functional analysis. A similarity matrix was generated based on the GO terms associated with the enriched modules. Each entry in the matrix represents the degree of functional similarity between modules in the two networks. Modules with highly similar biological functions between the two networks were identified. 25 modules were found to be enriched with the same GO terms and there are 52 modules with at least 50% similarity in GO terms.

Furthermore, to ascertain the representation of each biological process within the network, we quantified the number of modules associated with a particular biological process based on GO enrichment analysis. We determined the quantity of modules representing each biological process, offering a perspective on the distribution of these processes throughout the network (Figure 3).

4. Discussion

4.1. Gene Coexpression Network Alterations in Alzheimer’s Disease: Structural and Connectivity Insights

We used the ROSMAP cohort because it represents a broad-ranging spectrum of the older adult population compared to most autopsy series, which typically rely on selective recruitment criteria and individuals seeking medical care for their symptoms. This provides more representative and generalizable results. Additionally, we selected the DLPFC for this study due to its involvement in various higher cognitive functions, such as decision-making, working memory, and executive control, which are particularly relevant in AD, as impairments in these functions are hallmark features of the disease.

Previous transcriptomic analysis has revealed gene regulatory interactions occurring in brain tissues or regions of healthy and diseased individuals. For example, the impact of Alzheimer’s disease on different cortical regions has been investigated by gene co-expression network analysis and module and trait network analysis [43,44,66]. Here, two models of gene-to-gene correlations in the DLPFC of individuals with Alzheimer’s disease and a control group composed of older individuals were analyzed to examine alterations in coexpression patterns associated with Alzheimer’s pathology. Networks were constructed following restrictions that made them comparable. Both networks maintain the same number of connections and have a very similar size and diameter.

The K-S test indicated no significant statistical difference between the global degree distributions of the two transcriptional profiles. This suggests that the overall degree distribution remains largely unchanged. However, while the degree distributions appear statistically similar, differences may still emerge at mesoscopic and micro-level network properties. It is important to consider that both distributions correspond to the same anatomical regions in aged individuals, which may influence their structural and functional characteristics. Nonetheless, it is noteworthy that the AD network exhibits a higher number of highly co-expressed genes (see Figure 4). In contrast, the control network shows a more balanced distribution of connectivity where fewer nodes have a high degree of connections, indicative of a network that is less dependent on a small set of hubs for its overall cohesion.

Only around two-thirds of the gene-gene interactions are shared between the AD and the control network. Gene interaction patterns overlap substantially, but there is still about 31.61% difference in overall connectivity. However, this difference is key in differentiating between the normal and AD-related transcriptional regulation landscapes. The shared edges may represent preserved connections essential for basic brain function or transcriptional network neighborhoods less affected by AD pathology.

Additionally, the control network exhibits a slightly higher transitivity (0.6253) compared to the AD network (0.5956). While this difference is modest, the reduction in transitivity in the AD network may be showing a lower degree of local cohesion relative to the control network.

4.2. Functional Insights into Co-Expression Changes in the AD Network: Linking Epigenetic, Cytoskeletal, Immune, and Post-Transcriptional Pathways

In the context of co-expression networks, node degree can be interpreted as an indicator of a gene’s influence within the regulatory program. In the AD network, 13 hub genes were identified that are absent in the control network. These genes, listed in order of degree, include KRT6B, COL6A5, PDCD5P1, H3C9P, FCRL1, ONECUT3, CCL18, H2BC14, H2AC4, CROCC2, SLC25A24P1, PDE4DIP and FCRL3 (See Supplementary Table S3 for further details). A brief description of these genes and their relations with AD is listed in TableS3. Among these, KRT6B, located on chromosome 12, stands out as the most highly connected gene in the AD network, with a degree of 207.

From the 13 hub genes here, KRT6B, COL6A5 and CROCC2 are structural proteins and cytoskeletal components; H3C9P, H2BC14, H2AC4 and ONECUT3 are histones or genes related to chromatin regulation; FCRL1, FCRL3 and CCL18 are related to the immune system and inflammation; PDCD5P1 and PDE4DIP are pseudogenes related to neurodegeneration and cell death and SLC25A24P1 is a pseudogene of SLC25A24, involved in ATP transport and energy metabolism.

Interestingly, we can find 4 pseudogenes in the list of AD-specific-hub genes (SLC25A24P1, PDE4DIP, SLC25A24P1 and H3C9P). Studies have demonstrated that pseudogenes may participate in gene transcription post-transcriptional gene regulation and epigenetic regulation [67]. Other computationally predicted pseudogenes have been proposed to be involved in the etiology of AD and other NDs [68] and other pseudogenes have been identified as hubs in AD networks [43]. These findings suggest that pseudogenes may act as regulatory elements. Investigating the specific roles of these pseudogenes could provide deeper insights into their contributions to AD pathophysiology and highlight novel therapeutic targets for future research.

Other interesting gene is CCL18. It is known that chemokines and their receptors can promote AD pathology by inducing the production of A$\beta$ and some chemokines and receptors are also involved in tau phosphorylation. CCL18, a well described chemokine was found to have enhanced co-expression with other genes in the AD network. Given its enhanced co-expression with other genes in the AD network, CCL18 might act as a mediator or regulator of neuroinflammatory responses. For example, the chemokine-receptor interaction could stimulate pathways that promote the production or A$\beta$ accumulation [69].

We investigated the biological processes associated with the first neighbors of the hubs specific to the AD network. On average, changes in the degree of AD-specific hub genes affect 30.15 genes. Among these, CROCC exhibited the greatest increase in co-expression within the AD network, showing interactions with 66 newly co-expressed genes (See Figure 7a).In contrast, PDE4DIP displayed the fewest new co-expression connections, with only three additional genes.Furthermore, we identified a set of newly co-expressed genes that were consistently gained across all AD-specific hubs. Notably, 83 genes were found to be commonly co-expressed with hubs in the AD network. Functional enrichment analysis of these newly associated neighbors revealed four main functional groups See Figure 7b). The first group consists of processes related to nucleosome organization and assembly, primarily involving several histone genes. A key mechanism of epigenetic regulation in gene expression is the dynamic organization of chromatin structure, and epigenetic aberrations have been implicated in transcriptional alterations associated with AD [70].

The second group, which can be categorized as cytoskeletal organization and keratinization, is functionally related to the first group and includes genes from the KRT family as well as FLG. There is evidence suggesting that KRT family genes may play a role in neurological disorders, including dementia [71]. Research has demonstrated that the dysregulated expression of Keratin 9 is a consequence of AD pathology. However, its extensive interaction network implies additional, as yet unidentified downstream effects that may contribute to the initiation and progression of the disease [72]. Furthermore, keratin intermediate filaments have multiple roles in signaling pathways, inflammation, and various disease states [73]. Interestingly, the high-degree gene KRT6B has been implicated in the regulation of epithelial–mesenchymal transition (EMT) and immune-related pathways, including macrophage differentiation and polarization, as demonstrated in studies of bladder cancer [74].

This leads us to the third functional group, which includes antimicrobial and humoral immune response. In the past, research on inflammation-driven neurodegeneration and the immune role of A$\beta$ has given rise to the “Antimicrobial Protection Hypothesis” of AD. This hypothesis suggests that $\beta$-amyloid deposition represents an early innate immune response to actual or misinterpreted immunological threats. Initially, A$\beta$ captures and neutralizes invading pathogens by forming $\beta$-amyloid deposits. The subsequent fibrillization of A A$\beta$ activates neuroinflammatory pathways, aiding in pathogen clearance and $\beta$-amyloid removal. However, in AD, prolonged activation of this response results in persistent inflammation and neurodegeneration [75].

Finally, the fourth group is related to miRNA-mediated post-transcriptional gene silencing, a key pathway in the regulation of mRNA stability and translation, involved in fine gene control and adaptive response processes. The first three groups of metabolic pathways are mediated by at least one gene, and give rise to an explanation of the interpolation between at least two hypotheses of Alzheimer’s neuropathogenesis, including the antimicrobial protection hypothesis [70,76] and the neuroinflammation hypothesis [77]. The antimicrobial humoral response and the keratinization pathways of genes are connected by SPINK5, a serine protease inhibitor [78] and KRT6A, one of the 27 different type II keratins expressed in humans [79], while the immune response and nucleosome assembly pathways are at least connected by the H2BC10 histone.

4.3. High Betweenness Genes only Present in the Disease Network Are Involved in Diverse Biological Pathways

In the AD network, 44 genes high betweenness centrality genes were identified that were not central in the control network (See Supplementary Figure S4). Ten genes with higher betweenness were LOC124903002, RAB3A, SNAP91, BEX1, INPP5F, NECAP1, INA, ATP5F1B, NSG2 and SYP. LOC124903002, an uncharacterized gene, was the gene with higher betweenness. The gene appears in 5454.812 shortest paths between pairs of other nodes. Pathway enrichment analysis revealed that these AD-specific high-betweenness genes are primarily associated with processes related to synaptic function. Specifically, enriched pathways include the regulation of short-term neuronal synaptic plasticity, the synaptic vesicle cycle, synaptic vesicle maturation, the regulation of synaptic vesicle exocytosis and the establishment of vesicle localization. Beyond these, 20 other pathways were identified, all related to vesicle dynamics, exocytosis, transport, and regulatory mechanisms. It’s notable that high-betweenness centrality genes are enriched in synaptic processes, as synaptic processes are essential for sustaining the flow of information within the brain metabolism. These genes, acting as “bridges” or nodes of elevated information flow within the network, are playing a structurally analogous role to that of synaptic functions, which similarly manage and modulate information transmission across neural circuits.

In contrast, there are 39 genes that are serving as high information-passing nodes in the control network but are not shown to be highly mediated genes in the pathology network. Two nodes, KIFBP and ACTR2 are not even present in the AD network, due to their poor statistical dependency with other genes. These “missing” genes did not enrich any gene ontology as a group. The 10 higher betweenness centrality genes that are lost in the AD network as high degree-brokering genes include MAPK9, SYT3, PACSIN1, GNG3, FBXL19, AMPH, SERPINI1, RUSC2, AP2M1, and KIFBP. The Mapk9 mitogen-activated protein kinase 9 is the most highly mediated gene lost in the AD network. This gene has functions in the transduction of environmental signals to the transcriptional machinery, a function similar to those of other MAP kinases, and it phosphorylates the c-Jun protein [80], an essential component of the AP-1 transcription factor, highlights its involvement in stress-activated pathways, such as apoptosis, inflammation, and neuroprotection. The loss of MAPK9 (betweenness in control 5714.584, betweenness in AD network 738.7555) as a high betweenness gene in the AD network is particularly relevant because it implies a disruption in a key signaling pathway involved in cellular stress responses and transcriptional regulation, probably having extensive implications of network reconfiguration in AD and the potential functional vulnerabilities that may emerge from the disruption of signaling intermediates.

One of these genes, SERPINI1 (betweenness of 2179.33 and 570.30 in control and AD, respectively) encodes neuroserpin, an inhibitory serine protease crucial for regulating plasmin and other processes related to neuronal homeostasis. The aggregation of mutant neuroserpin polymers—resulting from the misfolding of neuroserpin—in the brain is associated with the formation of abnormal intracellular inclusions, leading to neuronal degeneration and dementia [81,82]. Nevertheless, there is evidence that neuroserpin is neuroprotective, as it reduces infarct volume and protects neurons from ischemia-induced apoptosis [83]. Neuroprotective effects may result from its role in inhibition of excitotoxicity, inflammation, and blood-brain barrier disruption that occur after acute ischemic stroke [84]. The node corresponding to SERPINI1 experienced a reduction of 1,609 in the number of shortest paths that pass through it, maybe reflecting a loss of neuroprotective capacity and contributing to the neuronal degeneration characteristic of AD.

Another gene that lost its status as a high-intermediate node was SYT3. Synaptotagmins are a class of Ca2+-dependent synaptic vesicle membrane proteins prevalent in the brain, potentially involved in membrane trafficking [85]. SYT3 is present in presynaptic terminals. It is expressed throughout the brain and is enriched in the brainstem, cerebellum and hippocampus [86]. Ischemic conditions have been reported to upregulate SYT3 expression [87]. Given the link between Alzheimer’s disease (AD) and impaired cerebral circulation—characterized by reduced oxygen and nutrient supply— Chum et al. [88] proposed that cerebrovascular miRNAs associated with AD regulate mRNA targets crucial for endothelial cell functions such as angiogenesis, vascular permeability, and blood flow regulation. Their findings revealed that several miRNAs associated with AD also regulate both cardiovascular and neuronal pathways, with SYT3 identified as a gene targeted by more than four miRNAs. This overlap shows a probable relationship between vascular and neuronal health in AD. Genes such PACSIN1, AMPH, RUSC2 and KIFBP participate in mechanisms associated with cytoskeleton dynamics and intracellular transport [89,90,91,92]. The absence of these genes as high-betweenness nodes indicates a disturbance in their function as mediators within the network, potentially hindering the coordination of cytoskeleton remodeling and intracellular trafficking pathways. Such abnormalities may lead to deficiencies in neuronal architecture, synapse functionality, and transport mechanisms, which are essential processes in neurodegenerative disorders such as AD.

4.4. Network Functional Analysis Show Gene Modular Rearrangement

The idea that sets of genes carrying out related tasks ought to be strongly correlated is supported by the clustering of genes in highly associated modules [93]. It has been seen that in gene co-expression models, different gene modules may correspond to distinct biological processes [94]. This observation aligns with the modular structure of the networks analyzed, where functionally related genes tend to co-cluster within specific modules. The AD network showed more modules than the control network. This statement should be taken with caution because it is known that Infomap tends to generate small size modules. Also, the modularity of the AD graph is higher (0.28) compared to the control graph (0.20), showing that the AD network has a slightly clearer community division. The control network has a slightly higher clustering coefficient (0.625) than the AD network (0.596). This indicates that, on average, nodes in the control network tend to form triangles (groups of three interconnected nodes) more frequently than in the AD network. The latter has more defined modules, but with less local cohesion, while the control graph has a less modular structure, but with a higher density of local connections. Under disease conditions, interactions between modules are more isolated, but connections within modules may not be as densely interrelated.

NMI indicates a 45% of similarity between the two graphs partitions. The similarity matrix based on the correspondence of shared genes between these modules revealed a reorganization of genes across various modules. While the modules in the AD and control networks are not directly equivalent, there is a degree of similarity between certain modules in both networks. However, a marked difference in the arrangement of gene membership across modules is apparent (see Figure 6a). In parallel, the biological similarity matrix for the communities displays a more conserved pattern. Although some modules show distinct biological functions, indicating possible divergence in their roles across different biological contexts, the functional organization remains relatively preserved despite the substantial reorganization of genes across modules (see Figure 6b).

Also, most biological processes exhibit notably higher module representation in the AD network compared to the control, where we can see that biological processes are distributed across a greater number of modules than in the control network (See Figure 3). This broader distribution may be the result of adaptation in the transcriptional program in response to the disease. Gene expression patterns are reorganized and compartmentalized to accommodate the altered cellular environment.

Axonogenesis and axon development display the greatest divergence in module representation, being distributed across 10 modules in the control network and 19 modules in the AD network, pointing out a possible emphasis on neuronal structural processes in AD. Other biological processes associated with AD, such as cognition and learning or memory, are also distributed across more modules in the AD network. Overall, there is a tendency in the AD network for biological processes to be represented in a greater number of modules. Notably, this phenomenon of spreading is also observed in processes such as rhythmic processes, axon guidance, learning and memory, response to amyloid-beta, and differentiation. Nevertheless, although fewer in number and of lesser magnitude, some processes exhibit the opposite trend. These include the regulation of protein secretion, positive regulation of secretion, amino acid and amide transport, as well as other transport- and movement-related molecular processes.

The reorganization of biological processes into a larger number of modules in the AD network may have several explanations, based both on AD biology and on how gene networks respond to stress and pathological changes. In AD, cells, especially neurons and glial cells, undergo significant changes in their environment due to factors such as the accumulation of beta-amyloid plaques, NFTs and neuroinflammation [95]. These changes generate a chronic cellular stress environment, which can trigger adaptive responses in gene expression. To cope with this adverse environment, the transcriptional program could be reorganized, promoting the activation and redistribution of genes co-expresssion. For instance, modules associated with processes such as axonogenesis, axon development, neuron projection guidance, and dendrite morphogenesis—critical for repairing and reorganizing damaged neuronal connections—show increased representation. This phenomenon reflects a possible compensatory or adaptive strategy to neurodegeneration.

In addition, greater representation of processes such as cognition and learning and memory in multiple modules in the AD network suggests that the disease affects entire sets of processes related to these high cognitive functions.

Differences in the transcriptional regulatory landscape between AD and healthy aged individuals are more subtle than in other diseases such as cancer, where we can see a clearer disruption in the modularity and coexpression [96]. The overall differences in topology point to changes at the mesoscopic level in the rearrangement of genes and molecular functions in co-expression modules in AD.

Moreover, because our correlations are based on the statistical dependence of mutual information, we can detect nonlinear relationships. However, these interactions lack directional inference, meaning we cannot determine whether they are associated with upregulation or downregulation. Additionally, AD is a disease with high comorbidity, a factor not explicitly accounted for in our analysis. More regions in the brain should be explored to obtain a more solid landscape of the complex transcriptional mechanism of the brain in AD, specially areas more affected by proteinopathy. Our findings should further be viewed as in-silico genomic-data-driven hypotheses that require additional experimental validation, given the statistical nature (gain/loss of correlation) of the differential connectivity relationships.