Identification of tryptophan metabolism-related genes in immunity and immunotherapy in Alzheimer’s disease

Raw data processing

GEO database was utilized using the following search criteria: Series: GSE5281, GSE37263, GSE106241, GSE132903, and GSE63060; Platform: GPL570, GPL5175, GPL6947, GPL24170, and GPL10558. The datasets GSE5281, GSE37263, GSE106241, and GSE132903 were utilized as training and test groups, while GSE63060 served as the independent test group. The search strategy employed was as follows: (‘Alzheimer’s disease’ [MeSH] mRNA [All Fields] and normal) AND (‘Homo sapiens’ [Organism] AND ‘Non-coding RNA profiling by array’ [Filter]). To identify relevant tryptophan metabolism-related genes, the Molecular Signatures Database (MSigDB) was employed, resulting in the retrieval of 40 genes (Appendix 1).

Analysis of DEGs

Transcription data obtained from Perl (https://github.com/Perl) were processed to obtain precise mRNA data. The obtained IDs were converted into corresponding gene names. Following this, data standardization was performed on GSE5281, GSE37263, GSE106241, and GSE132903 using the “normalize Between Arrays” function from the “limma” package. Principal Component Analysis (PCA) was conducted using the “factoextra” package to examine the variability in the data. To identify Differentially Expressed Genes (DEGs) between AD and non-demented controls (ND), statistical analysis was performed. The DEGs were selected based on the criteria of |Fold Change (FC)| > 1 and p-value < 0.05. To visually represent the significantly deregulated genes, a heatmap was generated using the “ggplot2” and “ComplexHeatmap” packages.

Immune cell infiltration

The immune cell components in adipose tissue were analyzed via CIBERSORT. We built barplot and corplot with the limma package to show the results of immune cells.

Cluster analysis

Cluster analysis was performed using the limma and ConsensusClusterPlus packages. By setting the clustering variable (k) to 2, we observed a strong correlation within each cluster and a weak correlation between clusters. The tryptophan metabolism-related genes associated with prognosis were classified into cluster 1 and cluster 2 based on this analysis. Additionally, a consensus score was calculated to further evaluate the results [12]. Furthermore, the limma package was employed to identify specific gene changes between subtypes and tissue types, allowing for the detection of differential gene expression patterns [13].

Enrichment analysis

To investigate the biological functions and pathways associated with the DEGs, we utilized Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). The DEGs were analyzed to determine their involvement in various biological processes, molecular functions, and cellular components using R. Specifically, we obtained the “c5.go.bp.v7.5.1.symbols” gene sets from MSigDB to perform the analysis.

To quantify the process scores, we employed the GSVA (Gene Set Variation Analysis) method, which is implemented in the “GSVA” package. This allowed us to compute the scores associated with different biological processes regulated by the DEGs. By utilizing this approach, we gained insights into the functional implications of the tryptophan metabolism-related genes in AD.

Co-expression gene identification

To explore the genetic processes involved in the development of AD, we employed the WGCNA algorithm. This algorithm allows for the clustering of genes into distinct modules and the examination of correlations between these modules and disease features. We utilized the “WGCNA” package to construct a co-expression network. To ensure robust results, we focused on genes with the highest 25% variance from the GSE132903 dataset. The dynamic cutting tree approach with a threshold of 0.25 was employed to merge modules. Several criteria were employed during the construction of the co-expression network. Firstly, the soft threshold power (β) was determined based on the scale-free topology requirement, aiming for a high independence value (R^2 = 0.85). The selection of the soft threshold was carried out using the select Soft Threshold function. Additionally, a minimum of 30 genes was set as the threshold for each module.

To identify potential relationships between the modules and clinical variables of the patients, Pearson correlation analysis was performed. This analysis allowed us to assess the degree of correlation between module expression patterns and clinical characteristics, providing valuable insights into the underlying mechanisms of AD.

Tryptophan metabolism-related genes identification

To identify the tryptophan metabolism-related genes, we performed an intersection analysis of the DEGs derived from major modules obtained through WGCNA, Gln, and cluster hubGenes. Visualization of the overlapping genes was achieved using Vnnmap. Furthermore, we investigated the biological processes and enrichment pathways associated with these genes. Separately, the datasets GSE5281, GSE37263, GSE106241, and GSE132903 were divided into training cohorts after identifying the hub tryptophan metabolism-related genes. To identify the hub DEGs, we utilized the “glmnet” package and selected the smallest lambda value as the optimal parameter. Subsequently, we calculated the DEGs’ predictive scores in each sample. To assess the diagnostic and discriminative potential of the tryptophan metabolism-related genes in AD and ND individuals, we conducted receiver operating characteristic curve analysis. The GSE63060 dataset was used for external validation.

Finally, we performed prognosis calculations on the test group by matching samples based on age-related clinical information. We also explored the correlation between these genes and age, providing valuable insights into their relationship with the aging process.

Drug-gene interactions

With the advancement of bioinformatics, the construction of biological models and the identification of efficient biomarkers has become more significant in the diagnosis and prevention of clinical disorders. Even if the biomarkers are established, the crucial issue is determining how to use them in the clinic. As a result, medication prediction based on successful indicators will be critical in the future prevention and treatment of AD. Validated biomarkers provide some reference for clinical treatment. Therefore effective drug prediction is very important. We used the DGIdb database (https://dgidb.genome.wustl.edu/) to make drug predictions for both the obtained hub genes and the intersection gene in the eXtreme Gradient Boosting (XGB) model.

Animals

Male APP/PS1 double transgenic mice (B6C3-Tg (APPswe, PSEN1dE9) 85Dbo/J) at the age of 12 months were obtained from Nanjing Junke Biotechnology Corporation, Ltd. (Nanjing, China). The experimental animal production license number was SCXK (SU) 2017-0003. The control group consisted of male C57BL/6J mice, also 12 months old, obtained from Hunan SJA Laboratory Animal Co., Ltd. (Changsha, China). The experimental animal production license number was SCXK (Xiang) 2019-0004. The mice were housed in a specific pathogen-free (SPF) animal room at Hunan University of Chinese Medicine. The animal reproduction and genotypic identification followed previously published methods [14].

Quantitative real time PCR (qPCR)

Rat hippocampal tissues were used for total RNA extraction using TRlzol (15596-026, Thermo Scientific, Shanghai, China). The extracted RNA was then subjected to reverse transcription using a reverse transcription cDNA kit (RR047A, Takara, Dalian, China) according to the provided instructions, resulting in cDNA synthesis. Two-step amplification of the target genes was conducted using SYBR Green (RR820L, Takara, Dalian, China) on a Real-time PCR instrument (T100, Bio-Rad, USA).

The amplification reaction procedure consisted of an initial denaturation step at 95° C for 5 minutes, followed by 40 cycles of denaturation at 95° C for 30 seconds and annealing/extension at 58° C for 30 seconds. Ct values, representing the cycle threshold, were recorded for each gene, and subsequent statistical analysis of the data was performed using the 2^-ΔΔCt method. The primer sequences used in the PCR amplification are presented in Table 1.

Western blotting

The hippocampus of mice was used for total protein extraction using RIPA lysate (cat. P0013B, Beyotime Biotechnology Co., Ltd., Shanghai, China) following the provided instructions. The protein concentration was determined using the BCA (23227, Thermo Scientific, Shanghai, China) method. Subsequently, the protein samples were loaded onto SDS-PAGE gels and subjected to electrophoretic separation. The separated proteins were then transferred to PVDF membranes (0.45 μm pore size, IPVH00010, Millipore Sigma Inc., Billerica, USA). To prevent non-specific binding, the membranes were blocked using 5% skim milk. The primary antibody (diluted at 1:1000) was incubated with the membranes overnight at 4° C, followed by incubation with the secondary antibody (diluted at 1:5000) at room temperature for 2 hours. Chemiluminescent signals were generated using an ECL chemiluminescence kit (32132, Thermo Scientific, Shanghai, China), and the blots were visualized using an imaging system (ChemiDoc™ XRS+, Bio-Rad, USA). The primary antibodies used in this study included Anti-PCCB Rabbit pAb (GB113237, Servicebio, Wuhan, China), Anti-Neurofascin Rabbit pAb (GB111351, Servicebio, Wuhan, China), and Anti-GAPDH mouse pAb (GB15002, Servicebio, Wuhan, China).

Immunofluorescence

Mouse brain tissue samples were fixed using 4% paraformaldehyde, followed by paraffin embedding. The paraffin-embedded tissue blocks were then sectioned into 4 μm thick slices using a paraffin microtome (HM355S, Thermo Scientific, Shanghai, China). To prepare the tissue for immunostaining, the sections underwent steps including deparaffinization, antigen retrieval, blocking, and incubation with primary and secondary antibodies, as described in previous studies [15]. Specifically, cell transparency techniques were employed, and antigen retrieval was performed to enhance the antigenicity of the tissue. After blocking to reduce non-specific binding, the sections were incubated with primary antibodies, including PCCB antibody at a dilution of 1:200 and Neurofascin antibody at a dilution of 1:200. Following the primary antibody incubation, the sections were incubated with appropriate secondary antibodies. Restaining and tablet sealing procedures were carried out according to established protocols. Full image scanning and subsequent statistical analysis were performed using the TissueFAXS imaging system (Tissue Gnostics GmbH, Austria), allowing for comprehensive evaluation of the stained tissue sections.

Statistical analysis

Statistical analysis was conducted using GraphPad Prism 8.0 Software. The data were expressed as mean ± standard deviation (SD). A two-group comparison was performed using the t-test. P-values less than 0.05 were considered statistically significant. For the identification of differentially expressed mRNAs, a threshold of FC greater than or equal to 2 and a p-value less than 0.05 were used for screening.

Data availability statement

The datasets generated during and/or analyzed during the current study are available in the Appendix.

DEG identification and principal component analysis

Among the 17 tryptophan metabolism-related genes examined, all showed significant differences between groups, except for ACMSD, GCDH, HADHA, HADH, MAOA, ALDH2, ALDH3A2, ALDH1B1, ALDH7A1, and ALDH9A1 (Figure 2A). These genes exhibited varying patterns of expression, with some genes clustering in the AD group and others in the control group. Specifically, the AD group showed increased expression of INMT, OGDH, MAOB, AADAT, IDO2, EHHADH, KYNU, ACAT1, TPH2, and CYP1A1. On the other hand, the control group exhibited increased expression of AOX1, CYP1A2, DDC, HAAO, OGDHL, TDO2, and ACAT2 (Figure 2B and Appendix 2).

Figure 2. Principal component analysis. (A) Tryptophan metabolism-related genes. (B) Expression of tryptophan metabolism-related genes in clusters.

Expression of tryptophan metabolism-related genes

To determine the chromosomal positions of the tryptophan metabolism-related genes, we performed calculations and visualized the data in a circular format (Figure 3A and Appendix 3). Additionally, in order to gain insights into the expression patterns of these genes, we conducted correlation analysis among them (Figure 3B, 3C).

Figure 3. Expression of tryptophan metabolism-related genes (TrpMgs). (A) Expression of TrpMgs on sequences. (B, C) The correlation between TrpMgs and related genes.

Immune cells

The immune environment plays a critical role in the initiation and advancement of AD. To assess the composition of immune cell populations in adipose tissue, we utilized CIBERSORT. Subsequently, we visualized the results using bar plots and correlation plots to depict the immune cell components (Figure 4A, 4B). Furthermore, to elucidate the relationship between the expression of the tryptophan metabolism-related genes and immune cells, we conducted correlation analysis (Figure 4C).

Figure 4. Expression of immune cells. (A, B) Expression of immune cells in different clusters. (C) Correlation between TrpMgs and immune cells.

Cluster analysis

By setting the clustering variable (k) to 2, we observed that the intragroup correlations were the strongest and the intergroup correlations were the smallest, suggesting that patients with AD can be classified into two distinct groups based on the expression of tryptophan metabolism-related genes (Figure 5A). Furthermore, we examined the expression patterns of these genes within the different clusters. Notably, TDO2, HAAO, ACAT2, and DDC did not exhibit significant differences between the two groups (Figure 5B, 5C). Additionally, based on the results of principal component analysis (PCA), patients with varying risk profiles were successfully segregated into two distinct groups (Figure 5D). Moreover, we investigated the immune cell infiltration patterns in relation to the different clusters, shedding light on the immune microenvironment associated with each group (Figure 5E, 5F).

Figure 5. Cluster analysis. (A) Consensus clustering matrix. (B, C) Expression of the TrpMgs in different clusters. (D) PCA. (E, F) Immune cell infiltration of different clusters.

Analysis of functional enrichments

Enrichment analysis was performed using the tryptophan metabolism-related genes. The molecular function (MF) analysis revealed their involvement in various functions, including transcription coregulator binding and histone arginine and methyltransferase activity. The biological processes (BP) associated with these genes encompassed diverse functions such as the regulation of cell migration, Aβ clearance, regulation of neurotransmitter transport, tryptophan catabolic process to kynurenine, serotonin biosynthetic process, and glycolytic process. Additionally, the cellular component (CC) analysis indicated their localization in the cytosol (Supplementary Figure 1A).

Pathway analysis further elucidated the functional implications of these genes. The enriched pathways included the chemokine signaling pathway, neuroactive ligand receptor interaction, p53 signaling pathway, apoptosis, tryptophan metabolism, and serotonergic synapse (Supplementary Figure 1B). These findings provide valuable insights into the molecular mechanisms and pathways associated with tryptophan metabolism in the context of Alzheimer’s disease.

Building a co-expression network and module detection

To ensure a scale-free topology in the co-expression network, a soft-thresholding power was employed to establish an appropriate approximation (Supplementary Figure 2A). Subsequently, the genes with the highest variance were organized and integrated into seven distinct co-expression modules (Supplementary Figure 2B). To explore the relationship between module eigengenes and clinical characteristics, Pearson’s correlation analysis was performed (Supplementary Figure 2C). Notably, the turquoise module exhibited strong connections with the “Group” attribute, which represents the classification of AD and non-AD samples, indicating its significant association with the disease (Supplementary Figure 2D and Appendix 4). These findings shed light on the potential relevance of the turquoise module in Alzheimer’s disease.

Clustering co-expression network construction and module detection

To establish a scale-free topology approximation for the network, a soft-thresholding power was applied, as depicted in Figure 6A. By clustering the genes based on their variance, co-expression modules were formed, as shown in Figure 6B. Pearson’s correlation analysis was then employed to investigate the relationship between module eigengenes and clinical characteristics, as illustrated in Figure 6C. Notably, the module exhibited strong connections with the “Group” attribute, representing the classification of AD and ND samples, indicating its significant association with the disease, as depicted in Figure 6D and (Appendix 5). The genes with the highest variance were further organized into eight co-expression modules, as presented in Figure 6E. Among these modules, the turquoise module demonstrated a high level of connectivity with the “Group” attribute, suggesting a strong association with AD and ND samples, as shown in Figure 6F.

Figure 6. Co-expression module construction. (A) Soft threshold power mean connection and scale-free fitting index analysis. (B) Dendrogram clustering. (C) Heatmap of correlations between module eigengenes and clinical characteristics. (D) Gene scatterplot in the turquoise module. (E) Clustering of dendrograms according to dynamic tree cutting, the genes were sorted into distinct modules using hierarchical clustering with a threshold of 0.25. Each color represents a separate module. (F) Gene scatterplot in the turquoise module.

Developing a model for least absolute shrinkage and operator selection

The intersection of DEGs, grey module genes (from WGCNA analysis), and tryptophan metabolism-related genes resulted in overlapping genes, as shown in Figure 7A and (Appendix 6). The boxplots illustrate the residual expression patterns of these genes in AD samples, highlighting the differences between AD and control groups, as presented in Figure 7B. Additionally, the proportions of the four different modes (DEGs, grey module genes, tryptophan metabolism-related genes, and overlapping genes) exhibit distinct variations, as depicted in Figure 7C. The expression of predictive values from the four models at different stages demonstrates noticeable differences, as shown in Figure 7D. Notably, the tryptophan metabolism-related genes exhibit a satisfactory diagnostic capacity in distinguishing AD from control samples, with area under the curve (AUC) values of RRF: 0.924, SVM: 0.937, XGB: 0.950, GLM: 0.641, as presented in Figure 7E. Among these models, the XGB model demonstrates the highest accuracy and stability (Appendix 7).

Figure 7. Cluster construction of co-expression modules. (A) Identification of tryptophan metabolism-related genes with a Venn diagram. (B, C) Residual expression patterns. (D) Model trend chart. (E) AUC of train group.

Model validation

In the external validation dataset GSE63060, the tryptophan metabolism-related genes exhibit an AUC of 0.924 (95% CI 0.876-0.963), indicating their potential as diagnostic markers for AD (Figure 8A). Additionally, the correlation analysis between the five hub genes and age revealed that PCCB and FARSB were negatively correlated with age, while NFASC and TEAD1 were positively correlated with age (Appendix 8). However, it is important to note that the P values for FARSB, NFASC, and TEAD1 were greater than 0.05, suggesting that these correlations may not be statistically significant (Figure 8B).

Figure 8. (A) AUC of test group. (B) Analysis of the relationship between hub genes and age.

Drug-gene interactions

According to the XGB model, one of the hub genes predicted three drugs, and one of these drugs is Dexamethasone (Table 2). Furthermore, we conducted predictions for all interacting genes associated with drug-gene relationships, and the results can be found in Appendix 9.

Validation of the 5-gene signature in APP/PS1 transgenic mice

The qPCR analysis revealed that the mRNA levels of PCCB and NFASC were down-regulated in the AD group compared to the control group, which is consistent with the previous analysis. However, the expression profile of FARSB did not align with the previous results. Furthermore, there were no significant differences in the expression levels of TEAD1 and EZR between the AD group and the wild-type group (Figure 9A–9E). The Western blotting results (Figure 9F and Supplementary Figure 3) demonstrated that the expression of PCCB and NFASC in the cerebral cortex of APP/PS1 transgenic mice was significantly decreased compared to the wild-type group. Additionally, immunofluorescence staining showed decreased expression of PCCB and NFASC in the cerebral cortex and hippocampus of APP/PS1 transgenic mice (Figure 9G–9H).

Figure 9. Validation of 5 gene markers in APP/PS1 mice. (A–E) qPCR to determine the mRNA expression of PCBB (A), FARSB (B), EZR (C), NFASC (D), and TEAD1 (E). (F) Western blotting to determine the protein expression levels of PCCB and NFASC. (G, H) Immunofluorescence to assess the expression of PCCB (G) and NFASC (H) in the hippocampus and cerebral cortex.