Integrative analysis of DNA methylation and gene expression identifies genes associated with biological aging in Alzheimer's disease.

Introduction: The acceleration of biological aging is a risk factor for Alzheimer's disease (AD). Here, we performed weighted gene co-expression network analysis (WGCNA) to identify modules and dysregulated genesinvolved in biological aging in AD.
Methods: We performed WGCNA to identify modules associated with biological clocks and hub genes of the module with the highest module significance. In addition, we performed differential expression analysis and association analysis with AD biomarkers.
Results: WGCNA identified five modules associated with biological clocks, with the module designated as "purple" showing the strongest association. Functional enrichment analysis revealed that the purple module was related to cell migration and death. Ten genes were identified as hub genes in purple modules, of which CX3CR1 was downregulated in AD and low levels of CX3CR1 expression were associated with AD biomarkers.
Conclusion: Network analysis identified genes associated with biological clocks, which suggests the genetic architecture underlying biological aging in AD. Highlights: Examine links between Alzheimer's disease (AD) peripheral transcriptome and biological aging changes.Weighted gene co-expression network analysis (WGCNA) found five modules related to biological aging.Among the hub genes of the module, CX3CR1 was downregulated in AD.The CX3CR1 expression level was associated with cognitive performance and brain atrophy.

INTRODUCTION

Alzheimer's disease (AD) is the most common neurodegenerative disease, characterized by accumulation of amyloid beta (Aβ) plaques and neurofibrillary tangles in the brain. , AD develops from multiple factors, such as genetics, lifestyle factors, and aging. The actual cause of AD is still an open question, but the greatest known risk factor for AD is aging, although chronological age is not a sufficient marker of individual health outcomes and does not reflect sensitivity to aging‐related diseases. Over the past decade there have been many studies to identify aging biomarkers that measure biological aging processes, also called biological age or biological clocks. , , , , The acceleration of biological age has been reported to associate with biomarkers of AD including the presence of neuritic plaques, amyloid load, and the decline of global cognitive functioning, episodic, and working memory. , , , , ,

Telomere length (TL) and epigenetic clocks are the most promising biological clocks and have been reported to be associated with the risk of AD. Telomeres are repeated sequences at the end of chromosomes that are shortened with each cell cycle. , Telomere shortening was associated with aging, age‐related diseases, including AD, and cognitive performance in AD. , , , , Epigenetic clocks are based on the changes in DNA methylation levels of multiple cytosine‐guanine dinucleotide (CpG) sites. Age‐related changes in methylation levels have been studied , , , and may cause various age‐related diseases including AD. , , , Moreover, epigenetic clocks have been linked to AD pathology, and the acceleration of epigenetic clocks has been associated with neuritic plaques and amyloid load in the prefrontal cortex and the decline of cognitive function, episodic memory, and working memory. , , ,

With the rapid development of high‐throughput screening techniques, gene expression analysis has become increasingly popular in studying the molecular mechanisms of human diseases. Previous studies have identified genes that are significantly dysregulated in AD. , , , In addition, individual genes do not work alone but interact with other genes. , Network analysis of gene expression profiles may offer the potential for identifying dysregulated genes that play important roles in complex diseases. Weighted gene co‐expression network analysis (WGCNA) has been developed to identify co‐expression modules of highly correlated genes, evaluate associations between modules and clinical outcomes, and identify potential disease‐associated gene sets. , , , We used WGCNA to identify modules by constructing a gene network using the correlation patterns between gene expression levels and evaluated the associations between modules and biological clocks.

To date, several biological clocks have been proposed. , , , Although previous studies reported the association between biological clocks and AD, transcriptomic factors influencing biological aging in AD are unknown. The objective of this study is to identify pathways and genes associated with biological aging and AD biomarkers in the AD continuum by integrating RNA profiling, DNA methylation, telomere lengths, cognition, and magnetic resonance imaging (MRI) scans from 551 non‐Hispanic White participants in the Alzheimer's Disease Neuroimaging Initiative (ADNI). First, we used DNA methylation data to estimate DNA methylation–based epigenetic clocks , and DNA methylation–based TL. In addition, we used TL measured using quantitative polymerase chain reaction (qPCR) to compare to DNA methylation–based estimated TL. After estimating biological clocks, we performed WGCNA on microarray messenger RNA (mRNA) expression profiles to identify modules (gene sets) as significantly associated with biological clocks. Then we identified hub genes in a module and performed differential expression analysis and association analysis of hub genes with AD diagnosis and AD biomarkers.

METHODS

Study participants

The data used in the study were downloaded from the ADNI database (http://adni.loni.usc.edu/). ADNI is a longitudinal study with more than 50 sites across the United States and Canada launched as a public‐private partnership in 2004, and details have been described previously. , In this study, a total of 551 ADNI participants including 181 cognitively normal older adults (CN), 272 mild cognitive impairment (MCI), and 93 AD were used for analysis.

Gene expression profiling

Gene expression profiling from blood samples was performed on the Affymetrix Human Genome U219 Array (www.affymetrix.com, Santa Clara, CA, USA). Raw microarray expression data were processed according to the standard quality control (QC) procedures described previously. In brief, all probe sets were annotated with the reference human genome (hg19), and raw expression values were normalized. A sex discrepancy check was performed using expression levels of sex‐specific genes. After quality control, the RNA expression profiles contained 21,150 probes. To adjust for the confounding effects, we used batch effects and RNA integrity number (RIN) values as covariates.

Telomere length measurement

Telomere length (or TL) measurement was performed on DNA from blood or a buffy coat sample of ADNI participants using a real‐time qPCR) assay. TL was measured at the same visits with DNA methylation profiling and quantified by comparing the amount of the telomere amplification product (T) to that of a single copy gene (S). The relative telomere to single copy gene (T/S) ratio was calculated to yield a value that correlates with the average TL, and the T/S ratio was adjusted for experimental variables to correct for batch effects. The TL acquisition and quality control have been described in detail previously.

RESEARCH IN CONTEXT

Systematic review: The authors reviewed the literature using a PubMed and Google Scholar search. There is increasing evidence that the acceleration of biological aging has been reported to be associated with Alzheimer's disease (AD). It is, therefore, possible that expression levels of biological clock–associated genes could be associated with cognition and AD biomarkers.

Interpretation: This is the first study to show that biological clocks were associated with gene co‐expression network modules, whose hub genes were also associated with cognition and brain atrophy, providing the genetic architecture underlying biological aging in AD.

Future directions: Further replication and functional studies in larger independent cohorts and animal models should be performed to investigate the mechanistic roles of CX3CR1 in biological clocks and AD pathology.

DNA methylation profiling

DNA methylation was profiled in blood or buffy coat samples using Illumina EPIC chips (Illumina, Inc., San Diego, CA, USA) according to the Illumina protocols. A detailed protocol describing DNA measurement has been described previously. In brief, genomic DNA samples were obtained from NCRAD (National Centralized Repository for Alzheimer's Disease and Related Dementias), and normalization and quality control were conducted. One sample was excluded because it had no CpG calls, and four samples were excluded because ≥1% of the CpG sites had a detection P‐value > .05. The remaining 1915 samples were normalized using wateRmelon. After quality control and normalization, the data set of comprised beta values was used for further analyses.

Prediction of DNA methylation–based epigenetic clocks and telomere length

DNA methylation–based epigenetic clocks (DNAmGrimAge and DNAmAgeSkinBlood) , and TL (DNAmTL) were predicted using the online biological clock calculator (http://dnamage.genetics.ucla.edu/) (see Figure S1 for the correlation of biological clocks). DNAmGrimAge and DNAmAgeSkinBlood are biological clocks (in units of years), and DNAmTL predicts leukocyte TL (LTL in units of kilobase [kb]) based on DNA methylation at multiple CpG sites. It is well known that the abundance of different cell types in blood changes with age, , and DNA methylation–based epigenetic clocks and DNAmTL are also associated with age‐related changes in blood cell composition. , , To adjust for the effect of blood cell composition, we estimated blood cell counts using DNA methylation profiles based on the Horvath method and the Houseman method, and we included the following seven blood cell counts in the multivariable regression to avoid multi‐collinearity between blood cell counts: B cells, CD8+ T cells, monocytes, plasma blasts, exhausted cytotoxic CD8+ T cells, naïve CD8+ T, and naïve CD4+ T.

Cognitive performance

Cognitive composite scores were downloaded from the ADNI database. Comprehensive neuropsychological evaluations were used to calculate composite measures for memory (MEM) and executive functioning (EF). , Briefly, composite scores for memory included the Rey Auditory Verbal Learning Test, Mini‐Mental State Examination, Logical Memory I‐II, and AD Assessment Scale‐Cognitive Test. Composite scores for EF included Digit Span Backwards, Category Fluency, WAIS‐R (Wechsler Adult Intelligence Scale‐Revised) Digit Symbol Substitution, Trail Making Tests A and B, and the Clock Drawing Test.

Magnetic resonance imaging processing

Three‐dimensional T1‐weighted MR images were downloaded from ADNI, and all MR images were processed with the FreeSurfer software (http://surfer.nmr.mgh.harvard.edu/). The detailed procedure has been described previously , Mean cortical thicknesses in regions of interest (ROIs) were estimated using Desikan‐Killiany‐Tourville ROIs. ,

Statistical analysis

Construction of gene co‐expression network and module detection

The gene co‐expression network was constructed using the WGCNA package in R. WGCNA was used to identify the modules of highly correlated genes that are related to external traits. First, similarity matrices of gene expression were calculated using Pearson's correlation Sij = |cor(i,j)|. Next, according to the theory of the network construction algorithm described in , the adjacency matrix (aij) was constructed using the power function β to construct a scale‐free network aij = power (Sij, β) = |Sij|β. To satisfy scale‐free topology, scale independence and mean connectivity were tested using a gradient method ranging from 1 to 20, and an appropriate power value was used. In this study, the power of β = 8 was selected to construct a scale‐free network. Then the adjacency matrix (aij) was transformed into a topological overlap matrix (TOM). Finally, module detection was conducted by performing average linkage hierarchical clustering of 1‐TOM dissimilarity matrix; detected modules are designated by color. , , ,

Identification of modules associated with biological clocks

To identify functional modules in the gene co‐expression network, the module eigengene (ME) was summarized by the first principal component of a given module, and module‐trait relationships were estimated using the correlation between MEs and biological clocks. The associations between MEs and biological clocks were estimated after adjusting for age and sex. To adjust for multiple testing, we used false discovery rate (FDR) corrected P‐values. To evaluate correlation strength, we calculate module significance (MS), which is defined by the average of the absolute gene significance for all genes in a given module, where the gene significance is defined by the correlation between the trait and gene expression levels of the gene. The relationships between modules were estimated by correlations between MEs, and hierarchical clustering of MEs was performed. In the WGCNA, modules with the highest MS are usually defined as the key module and selected for further analysis. , ,

Functional enrichment analysis of a module

To understand the biological relationships of genes within a module, gene‐set enrichment analysis (GSEA) was performed to identify enriched biological pathways. We selected genes in a module of WGCNA and performed GSEA with pathways from the Gene Ontology (GO) database (http://geneontology.org/). Pathways with FDR‐corrected P‐value < .05 were considered as significant.

Construction of a protein‐protein interaction (PPI) network and identification of hub genes in a module

We constructed a protein‐protein interaction (PPI) network of genes in a module. To identify experimentally validated interactions in the module, PPIs from the STRING database (https://string‐db.org/) were used. Interactions with a confidence score >0.7 were used for network construction. To identify hub genes of the module, we used the maximal clique centrality (MCC) algorithm. The MCC of each protein node was calculated using cytoHubba in the Cytoscape platform (https://cytoscape.org/). In this study, genes with the top 10 MCC scores were considered hub genes of the module.

Associations of expression levels of hub genes with AD and AD endophenotypes

First, we performed differential expression analysis of probes belonging to hub genes between diagnostic groups using analysis of covariance (ANCOVA) with age and sex as covariates. Then we performed a post hoc analysis of probes significantly associated with the diagnostic groups using AD endophenotypes. We used composite scores for memory and EF and entorhinal cortical thickness measured from MRI scans as AD endophenotypes. We performed linear regression analysis to investigate the association between expression levels of probes and AD endophenotypes. Age, sex, education, and intracranial volume (ICV) were used as covariates, if appropriate.

RESULTS

A total of 551 ADNI non‐Hispanic White participants were used for gene co‐expression network analysis (181 CN, 272 MCI, and 93 AD) and demographic information is presented in Table 1. Age and sex were significantly different between diagnostic groups, whereas education and RIN values did not show any significant differences between diagnostic groups. As expected, the number of apolipoprotein E (APOE) ε4 alleles was greater in MCI and AD.

TABLE 1

Demographic characteristics of the study participants

	CN	MCI	AD	P
N (N = 551)	186	272	93
Age	76.2 (6.44)	72.6 (7.72)	77.1 (7.91)	4.66 × 10–9
Sex (female/male)	92/94	120/152	35/58	1.64 × 10–1
Education	16.4 (2.71)	16.2 (2.71)	15.9 (3.03)	4.43 × 10–1
Number of APOE ε4 alleles (0/1/2)	137/45/4	157/95/20	30/46/17	1.52 × 10–10
RIN	6.90 (0.52)	7.01 (0.52)	6.98 (0.64)	1.35 × 10–1

Abbreviations: AD, Alzheimer's disease; APOE, apolipoprotein E gene; CN, cognitive normal; MCI, mild cognitive impairment; RIN, RNA integrity number.

Note: Boldface indicates P‐values less than .05. The mean and standard deviation of age, education, RIN values are shown in the table.

Construction of gene co‐expression networks and identification of modules

To identify gene modules, the gene co‐expression networks were constructed using WGCNA. To ensure a scale‐free network, scale independence and mean connectivity were tested using the gradient method, with the power value ranging from 1 to 20, and the power of β = 8 was selected. After the construction of the scale‐free network, the modules were detected using a hierarchical clustering dendrogram (Figure S2). As the result, 13,302 probes were grouped into 21 modules using the average linkage hierarchical clustering, and 7848 un‐clustered probes (included in the gray module) were removed for subsequent analysis. Twenty‐one modules were labeled by color, and the number of genes in each module is shown in Table S1.

Identification of gene co‐expression modules associated with biological clocks

The module‐trait relationships of 21 modules were tested by calculating the correlation between the MEs and biological clocks. Among the 21 modules, five modules (purple, blue, brown, tan, and magenta) were associated significantly with TL, DNAmTL, and DNAmAgeSkinBlood after adjusting for multiple testing using FDR correction (Figure 1A). In particular, the purple module was significantly associated with TL (r = −0.19, FDR corrected P = 1.59 × 10–4), DNAmTL (r = −0.22, FDR corrected P = 6.25 × 10–6), and DNAmAgeSkinBlood (r = 0.2, FDR corrected P = 6.23 × 10–5) (see Tables S2, S3, and S4 for module‐trait relationships in each diagnostic group). MS is defined by the average of the absolute gene significance for all genes in a given module, where the gene significance is defined by the correlation between the trait and expression levels of the gene. WGCNA usually defined the module with the highest MS as the key module, and it was selected for further analysis. , , In this study, the purple module has the highest mean MS across the biological clocks (Figure 1B). In addition, the relationships between modules were tested. As shown in Figure S3, modules were divided primarily into two clusters according to the correlations between MEs. The results showed that, of the five modules that are significantly associated with biological clocks, the tan and magenta modules were clustered with the purple module, whereas the brown and blue modules were not. Therefore, based on the relationships between modules, further analysis was performed on the purple module, a key module, as well as two modules, blue and brown, which did not cluster with the purple module.

FIGURE 1

Identification of modules associated with biological clocks. (A) Relationship between modules and biological clocks. Heatmap displays the correlation coefficient (r) between module eigengenes (MEs) and biological clocks and false discovery rate (FDR)–corrected P‐values. (B) Average of the module significance (MS) of biological clocks.

Gene set enrichment analysis of genes within the modules

To understand the biological relationship of genes in the modules, GSEA was performed with probes belonging to each module. The GSEA analysis with pathways from GO showed 50 enriched pathways for the purple module (Table 2). For the measured TL, the purple module was significantly enriched in a nuclear protein‐containing complex (FDR corrected P = 1.90 × 10–2). For the DNAmTL, the purple module was enriched in 11 pathways including leukocyte migration. Finally, for the DNAmAgeSkinBlood, the purple module was significantly enriched in 39 pathways including cell and intracellular transport. The GSEA results for the blue and brown modules are shown in Tables S6 and S7. In the blue module, 25 pathways were enriched for TL and DNAmAgeSkinBlood, including transferase complex, leukocyte mediate immunity, and cell activation (Table S6). In the brown module, five pathways were enriched for DNAmGrimAge and DNAmAgeSkinBlood (Table S7).

TABLE 2

Gene‐set enrichment analysis of the purple module for biological clocks

Biological age	Gene Ontology	N genes	NES	FDR P
TL	Nuclear protein–containing complex	6	−1.130	1.90 × 10–2
DNAmTL	Regulation of response to biotic stimulus	13	1.210	4.10 × 10–2
	Response to bacterium	8	1.200	3.60 × 10–2
	Granulocyte migration	6	−1.138	4.66 × 10–2
	Leukocyte migration	14	−1.141	4.67 × 10–2
	Positive regulation of intracellular signal transduction	16	−1.127	4.68 × 10–2
	Positive regulation of MAPK cascade	7	−1.132	4.75 × 10–2
	Myeloid cell differentiation	6	−1.128	4.80 × 10–2
	Neutrophil migration	5	−1.142	4.82 × 10–2
	Regulation of leukocyte migration	7	−1.132	4.83 × 10–2
	Protein dimerization activity	6	−1.119	4.58 × 10–2
	Identical protein binding	21	−1.101	4.72 × 10–2
DNAmAgeSkinBlood	Positive regulation of intracellular transport	6	2.405	2.95 × 10–4
	Regulation of intracellular transport	6	2.405	3.69 × 10–4
	Positive regulation of locomotion	10	2.422	4.92 × 10–4
	Positive regulation of phosphorus metabolic process	13	2.441	7.38 × 10–4
	Taxis	14	2.326	9.98 × 10–4
	Regulation of cellular component movement	17	2.456	1.48 × 10–3
	Locomotion	31	2.249	1.72 × 10–3
	Positive regulation of intracellular protein transport	5	2.258	1.94 × 10–3
	Regulation of intracellular protein transport	5	2.258	2.22 × 10–3
	Positive regulation of protein localization to membrane	5	2.135	3.99 × 10–3
	Regulation of protein localization to membrane	5	2.135	4.35 × 10–3
	Regulation of anatomical structure morphogenesis	7	2.104	4.36 × 10–3
	Regulation of vasculature development	5	2.090	4.47 × 10–3
	Positive regulation of protein phosphorylation	12	2.144	4.52 × 10–3
	Second messenger mediated signaling	9	2.047	5.50 × 10–3
	Calcium mediated signaling	9	2.047	5.87 × 10–3
	Cell migration	27	1.992	6.83 × 10–3
	Protein kinase b signaling	8	2.003	6.87 × 10–3
	Gliogenesis	6	1.981	7.41 × 10–3
	Regulation of intracellular signal transduction	23	1.928	1.12 × 10–2
	Positive regulation of protein kinase b signaling	5	1.902	1.33 × 10–2
	Intracellular transport	13	1.878	1.64 × 10–2
	Positive regulation of catalytic activity	14	1.869	1.72 × 10–2
	Organophosphate metabolic process	6	1.864	1.72 × 10–2
	Blood vessel morphogenesis	10	1.857	1.79 × 10–2
	Cell chemotaxis	8	1.841	1.90 × 10–2
	Regulation of cell death	13	1.798	2.45 × 10–2
	G protein–coupled receptor signaling pathway	14	1.798	2.52 × 10–2
	Negative regulation of cell death	8	1.800	2.57 × 10–2
	Positive regulation of chemotaxis	5	1.773	2.70 × 10–2
	Regulation of chemotaxis	5	1.773	2.79 × 10–2
	Growth	8	1.775	2.86 × 10–2
	Anatomical structure formation involved in morphogenesis	13	1.743	3.10 × 10–2
	Positive regulation of protein modification process	15	1.736	3.14 × 10–2
	Cellular component morphogenesis	8	1.745	3.16 × 10–2
	Biological process involved in symbiotic interaction	17	1.701	4.22 × 10–2
	Positive regulation of intracellular signal transduction	16	1.685	4.63 × 10–2
	Vasculature development	11	1.669	4.99 × 10–2
	Cell morphogenesis	8	1.671	5.00 × 10–2

Abbreviations: FDR P, false discovery rate (FDR)–corrected P‐value; N genes, the number of genes in gene sets; NES, normalized enrichment score.

Identification of hub genes in the modules

The interaction networks in each module were constructed using the PPI network. The hub genes within a network were determined by the MCC algorithm. The PPI network with the genes in the purple module consisted of 49 nodes and 114 edges (Figure 2A). Genes with the top 10 higher MCC scores were selected as hub genes of the purple module (Figure 2B and Table S5): CCR5 (C‐C motif chemokine receptor 5), CCL4 (C‐C motif chemokine ligand 4), CXCR6 (C‐X‐C motif chemokine receptor 6), CX3CR1 (C‐X3‐C motif chemokine receptor 1), NMUR1 (neuromedin U receptor 1), S1PR5 (sphingosine‐1‐phosphate receptor 5), CCL5 (C‐C motif chemokine ligand 5), PRF1 (perforin 1), IL2RB (interleukin 2 receptor subunit beta), and GZMB (granzyme B). The PPI network with the genes in the blue module consisted of 1439 nodes and 5752 edges, and the network with the genes in the brown module consisted of 1398 nodes and 7874 edges. The hub genes in the blue and brown modules and their first connected genes are shown in Figures S4 and S5.

FIGURE 2

Visualization of the protein‐protein interaction (PPI) network and hub genes. (A) PPI network with genes in the purple module. The blue nodes represent the genes and edges present the associations between genes. (B) 10 Hub genes with highest maximal clique centrality (MCC) scores and their first connected genes. The MCC scores are displayed in colors from red to yellow on the node.

Association of blood expression levels of hub genes with AD diagnosis and AD biomarkers

We performed differential expression analysis of probes belonging to the 10 hub genes in each module. In total, 14 probes in the purple module, 17 probes in blue, and 21 probes in the brown module were used for the differential expression analysis. Among the 14 probes in the purple module, one probe annotated on CX3CR1 gene showed a significant difference of expression levels between diagnostic groups (Table 3 and Figure S6). CX3CR1 was significantly downregulated in AD compared to MCI and CN (P = 3.27 × 10–2, AD‐MCI P = 1.62 × 10–3, AD‐CN P = 1.40 × 10–3, MCI‐CN P = 4.81 × 10–1). There were no differentially expressed genes between diagnostic groups in the blue and brown modules (Tables S8 and S9). We evaluated the associations with AD biomarkers for the differentially expressed gene in the purple module. The expression levels of CX3CR1 were significantly associated with composite scores for MEM (β = 0.301, P = 2.27 × 10–2) and for EF (β = 0.407, P = 5.22 × 10–3) (Figure 3A,B). Higher expression levels of CX3CR1 were associated with higher composite scores for both MEM and EF. In addition, higher expression levels of CX3CR1 were also significantly associated with thicker entorhinal cortical thickness (β = 0.153, P = 2.40 × 10–2) (Figure 3C).

TABLE 3

Differentially expressed probes and related hub genes in the purple module

probe ID	Gene	CN	MCI	AD	P	FDR P
11723048_at	CX3CR1	0.032	0.005	−0.081	2.33 × 10–3	3.27 × 10–2
11752490_x_at	NMUR1	0.000	0.013	−0.044	5.54 × 10–1	7.65 × 10–1
11740639_at	NMUR1	−0.013	0.018	−0.033	6.52 × 10–1	7.65 × 10–1
11729977_a_at	CXCR6	0.033	−0.017	−0.012	6.55 × 10–1	7.65 × 10–1
11747295_a_at	PRF1	0.006	0.010	−0.046	6.02 × 10–1	7.65 × 10–1
11730909_s_at	CCR5	0.043	−0.014	−0.040	2.74 × 10–1	7.65 × 10–1
11718983_x_at	CCL4	0.027	−0.005	−0.038	3.49 × 10–1	7.65 × 10–1
11752664_a_at	S1PR5	0.004	0.007	−0.031	8.35 × 10–1	8.74 × 10–1
11722635_at	IL2RB	0.006	0.018	−0.073	1.15 × 10–1	7.65 × 10–1
11732275_at	CCL5	0.018	0.004	−0.050	4.79 × 10–1	7.65 × 10–1
11718982_s_at	CCL4	0.021	−0.005	−0.025	3.84 × 10–1	7.65 × 10–1
11753810_a_at	CCL5	0.021	0.004	−0.057	3.13 × 10–1	7.65 × 10–1
11732276_x_at	CCL5	0.022	0.007	−0.069	1.99 × 10–1	7.65 × 10–1
11724900_a_at	GZMB	0.001	0.005	−0.018	8.74 × 10–1	8.74 × 10–1

Abbreviations: AD, Alzheimer's disease; CN, cognitively normal old adults; FDR P, false discovery rate (FDR)−corrected P‐value; MCI, mild cognitive impairment; P, P‐value.

Note: Boldface indicates FDR P‐values less than 0.05.

FIGURE 3

Associations between expression levels of CX3CR1 gene and Alzheimer's disease (AD) biomarkers. Expression levels of CX3CR1 versus composite scores for memory (A), executive functioning (B), and entorhinal cortical thickness (C).

DISCUSSION

Accumulating evidence suggests that acceleration of biological aging may contribute to increasing risk for AD. , , , , Here, we conducted WGCNA of transcriptomics profiles to identify modules and hub genes of the modules associated with biological aging in the continuum of AD (MCI, AD) using four types of biological clocks including measured TL, two DNA methylation–based clocks (DNAmGrimAge, DNAmAgeSkinBlood), and DNA methylation–based estimated telomere length (DNAmTL).

Gene co‐expression network analysis identified five modules significantly associated with biological aging, and among them the purple module showed the strongest associations with four biological clocks. GSEA of the purple module identified 50 significantly enriched pathways related to the acceleration of biological clocks. Our analysis revealed the enriched pathways related to leukocyte migration and intracellular transport. Leukocyte migration plays an important role in inflammation, which is linked to AD development, and leukocyte filtration in AD brain is stimulated by Aβ. The disruption of the intracellular transport of enzymes involved in the production of Aβ has been associated with the onset of AD. Furthermore, pathway analysis also revealed pathways related to cell migration, cell death, chemotaxis, and gliogenesis. Cell death and migration are related to Aβ, a key molecule in the pathogenesis of AD. The migration of monocytes into the AD brain begins with Aβ deposition and contributes to neuronal cell death that causes brain atrophy. Microglial chemotaxis is related to the concentration of soluble Aβ, and migration of microglia induced by Aβ constitutes a chemotactic response of neuroimmune brain cells. Gliogenesis in the adult brain regulates brain function and is involved in neurological conditions such as AD. ,

Therefore, genes in the purple module are likely linked to biological processes relevant to AD. Furthermore, according to the MCC algorithm, the top 10 genes with higher MCC scores were considered as hub genes of the purple module, and their expression was found to be downregulated in patients with AD compared with cognitively normal older adults (or CN). Among them, CX3CR1 expression was significantly downregulated in AD and was also significantly associated with cognitive performance and brain atrophy.

The CX3CR1 (CX3C chemokine receptor 1) gene is a protein‐coding gene that binds to CX3CL1 (CX3C chemokine fractalkine), also known as fractalkine. The chemokine is involved in the adhesion and migration of leukocytes and chemokine, and chemokine receptors play a role in the development of AD. Genetic mutations in the CX3CR1 gene are associated with the number of cells in lymphocytes, , leukocytes, , , and myeloid white cells. CX3CL1/CX3CR1, one type of chemokine and chemokine receptors, has been reported to play a role in AD. , In the central nervous system (CNS), CX3CL1 is highly expressed in neurons, and CX3CR1 is mainly expressed in microglia. In our study, the CX3CR1 gene was downregulated in AD, which is consistent with the previous study. Previous studies suggest that the interruption of the CX3CL1/CX3CR1 signaling, which plays an important role in neuro‐microglia interaction, has beneficial as well as detrimental effects on AD pathology. CX3CL1/CX3CR1 signaling has a neuroprotective effect in the CNS by reducing neurotoxicity and microglial activation , and has a beneficial effect on Aβ clearance. , , , , Conversely, CX3CL1/CX3CR1 signaling has also been reported to have a negative effect on AD pathology. Disrupting the CX3CL1/CX3CR1 signaling pathway by ablation of the CX3CR1 resulted in enhancement of tau accumulation. , , , , In addition, the disruption of CX3CL1/CX3CR1 signaling has been associated with the dysregulation of microglial responses, neuronal damage, , alleviation of Aβ neurotoxicity, and memory deficiency. ,

Among the hub genes in the purple module, five genes (CX3CR1, CCL4, CCR5, CCL5, CXCR6) are chemokine or chemokine receptors. Chemokines and their receptors have been reported to be related to neuroinflammation in AD. CCR5 (C‐C chemokine receptor type 5) is involved in neuroinflammation processes associated with AD , , and related to long‐term and spatial memory functions in mice. In addition, the absence of CCR5 modulates the expression of CCR2, which is related to the activation of astrocytes, which may lead to Aβ deposition. CCL4 (C‐C motif chemokine ligand 4), also known as macrophage inflammatory protein‐1β (MIP‐1β), binds to CCR5 and CCR8 (C‐C chemokine receptor 8). The expression of CCL4 is modulated by Aβ deposition.

This study is limited due to its modest sample size. Replication studies in independent cohorts with large sample sizes are needed to validate our findings. In addition, we used only non‐Hispanic Whites from the ADNI cohort to avoid confounding by population stratification. Further studies in racially diverse samples will be important for generalization of these findings to community‐based and ancestrally diverse populations.

In conclusion, we performed a co‐expression network analysis of blood gene expression profiles to identify modules and dysregulated genes associated with biological clocks in AD. We identified five modules related to biological clocks, where the purple module showed the most significant association with biological clocks. Among hub genes of the purple module, CX3CR1 was downregulated in AD and associated with cognitive performance and brain atrophy. The biological role of CX3CR1 in biological clocks merits further investigation. This study may contribute to the understanding of the genetic architecture underlying biological aging and its role in relation to AD pathophysiology.

CONFLICT OF INTEREST

Qingqin S. Li is an employee of Janssen Research & Development, LLC, an operating company of Johnson & Johnson, and holds company equity. The authors declare that they have no non‐financial competing interests. Author disclosures are available in the supporting information.

Supporting Information

Click here for additional data file.

Supporting Information

Click here for additional data file.