Abstract
Background
Alzheimer’s disease (AD) is the most common form of age-related
neurodegenerative disease. Unfortunately, due to the complexity of
pathological types and clinical heterogeneity of AD, there is a lack of
satisfactory treatment for AD. Previous studies have shown that
microRNAs and transcription factors can modulate genes associated with
AD, but the underlying pathophysiology remains unclear.
Methods
The datasets [35]GSE1297 and [36]GSE5281 were downloaded from the gene
expression omnibus (GEO) database and analyzed to obtain the
differentially expressed genes (DEGs) through the “R” language “limma”
package. The [37]GSE1297 dataset was analyzed by weighted correlation
network analysis (WGCNA), and the key gene modules were selected. Next,
gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG)
pathway enrichment analysis for the key gene modules were performed.
Then, the protein-protein interaction (PPI) network was constructed and
the hub genes were identified using the STRING database and Cytoscape
software. Finally, for the [38]GSE150693 dataset, the “R” package
“survivation” was used to integrate the data of survival time, AD
transformation status and 35 characteristics, and the key microRNAs
(miRNAs) were selected by Cox method. We also performed regression
analysis using least absolute shrinkage and selection operator
(Lasso)-Cox to construct and validate prognostic features associated
with the four key genes using different databases. We also tried to
find drugs targeting key genes through DrugBank database.
Results
GO and KEGG enrichment analysis showed that DEGs were mainly enriched
in pathways regulating chemical synaptic transmission, glutamatergic
synapses and Huntington’s disease. In addition, 10 hub genes were
selected from the PPI network by using the algorithm Between
Centrality. Then, four core genes (TBP, CDK7, GRM5, and GRIA1) were
selected by correlation with clinical information, and the established
model had very good prognosis in different databases. Finally,
hsa-miR-425-5p and hsa-miR-186-5p were determined by COX regression, AD
transformation status and aberrant miRNAs.
Conclusion
In conclusion, we tried to construct a network in which miRNAs and
transcription factors jointly regulate pathogenic genes, and described
the process that abnormal miRNAs and abnormal transcription factors TBP
and CDK7 jointly regulate the transcription of AD central genes GRM5
and GRIA1. The insights gained from this study offer the potential AD
biomarkers, which may be of assistance to the diagnose and therapy of
AD.
Keywords: Alzheimer’s disease (AD), differentially expressed genes
(DEGs), COX regression, transcription factors, miRNAs
Introduction
Alzheimer’s disease (AD) is the most common and complex neurological
disease in the world. Neurological conditions such as AD are the
leading cause of 70% of dementia worldwide ([39]Mayeux and Stern,
2012). Global Alzheimer’s cases are increasing year by year and are
expected to reach 78 million by 2030. While the pathogenesis of the
disease is not yet known, many believe the accumulation of
over-phosphorylated tau and Aβ plaques in neurofibrillary tangles are
the main causes of AD ([40]Andreasen et al., 1999; [41]Shao et al.,
2011; [42]Arvanitakis and Bennett, 2019). When extracellular amyloid
plaques build up in specific areas of the brain, they can lead to
amyloid vascular disease or neurodegenerative diseases ([43]Thanvi and
Robinson, 2006; [44]Jäkel et al., 2022). The neurofibrillary tangles
(NFTs) are huge paired intracellular helical strands of
hyperphosphorylated tau proteins, which induce neuronal and synaptic
loss ([45]Moloney et al., 2021). The predominant regions of the
pathological process underlying AD in human brain are the association
areas of the cerebral cortex and the hippocampus ([46]Wang et al.,
2010). Despite these extensive findings on both, there are few
effective drugs to improve and treat AD. Worse, many patients have to
wait a long time to be diagnosed with the disease. During this period,
the Aβ burden is significant and memory loss has already occurred
([47]Arvanitakis et al., 2019; [48]Long and Holtzman, 2019). Typical AD
has gone through a gradual and hidden development process, and there is
no specific detection method, which cannot be accurately diagnosed,
difficult to cure, difficult to control, and poor prognosis. Therefore,
it is urgent to explore new potential biomarkers for early diagnosis
and effective treatment of AD.
Alzheimer’s disease is the most common neurodegenerative disorder with
limited therapeutics, and AD is characterized by the formation of
plaques made by protein aggregates. Mounting studies have suggested
that targeting transcription factors holds promise for treating
neurodegenerative disorders including AD (Ping [49]Yang et al., 2022).
CMV promotor-driven transcription factor EB (TFEB), for example, is
injected with adeno-associated viral particles in targeted mice, which
are mainly localized in neuronal nuclei and upregulated lysosomes. This
resulted in decreased steady state levels of Aβ in APP proteins and
cerebral stromal liquid ([50]Xiao et al., 2015; [51]Song et al., 2020).
A growing body of evidence also suggests that the accumulation of
misfolding proteins in AD is caused by damage to macromolecular
autophagy and autophagy lysosomal pathways, making TFEB, which
regulates autophagic lysosomal pathway, a promising target for AD
treatment ([52]Zheng et al., 2021). We also found decreased expression
of the transcription factor Nrf2 and its NQO1, HO-1, and GCLC driver
genes, and changes in related Nrf2 pathways in AD brains. Nrf2
activation may provide cellular protection and prevent an increasing
number of diseases including neurodegenerative diseases. These avenues
of evidence point to the activation of the Nrf2 transcription factor as
a potential new therapeutic scheme for AD ([53]Osama et al., 2020).
Additional research shows that ATF6 transcription factor activation can
attenuate amyloidosis through BACE1 downregulation, and concomitantly
ATF6 overexpression significantly reduces Aβ1-42. The results of this
study suggest ATF6 may become a potential focus for targeted therapy of
AD. Interestingly, TBP and CDk7 also belong to the basal transcription
factor family like ATF6 among the 10 key genes selected in our study
([54]Geng et al., 2016; [55]Du et al., 2020).
Recently, progress has begun to clarify the physiological and
pathological roles of non-coding RNAs (ncRNAs) in a variety of
diseases, including cancer. ncRNAs do not have protein coding
functions, but are important regulators of many cellular processes
([56]Esteller, 2011). Of these, miRNAs are the most studied and have
become key players in the pathogenesis of AD involved in regulating key
growth regulatory pathways ([57]Swarbrick et al., 2019; [58]Toden et
al., 2021). miRNAs are small endogenous non-coding RNA molecules that
can inhibit or silence the expression of post-transcription genes. Many
miRNAs play an important role in post-transcriptional regulation and
are highly conserved ([59]Sand et al., 2009). Previous research has
found that numerous miRNAs, including miRNA-126a-3p, MicroRNA-455-3p,
miR-501-3p, and miRNA-101a-3p, have a role in the pathogenesis of AD,
showing that miRNAs are strongly linked to the development of AD
([60]Hara et al., 2017; [61]Kumar et al., 2017; [62]Kumar and Reddy,
2019; [63]Lin et al., 2022; [64]Xue et al., 2022). However, the
specific role and potential mechanisms of miRNA in AD pathogenesis are
unclear.
We here established a miRNA-mRNA, TF-mRNAs regulatory network by
integrating relevant TFs (transcription factors), miRNAs and mRNAs to
try to gain further insight into the potential functions of ncRNAs as
well as transcription factors in AD. This study may provide a better
understanding of the underlying pathogenesis of AD and potentially new
biomarkers for diagnosis and treatment of AD.
Materials and methods
Gene expression profile data collection
The gene expression datasets used in this investigation were collected
from the gene expression omnibus (GEO) database.^[65]1 This database
generated a total of 340 datasets on human AD. After a careful review,
four gene expression profiles ([66]GSE5281, [67]GSE1297, [68]GSE122063,
and [69]GSE150693) were selected. Among them, [70]GSE5281 was based on
the [71]GPL570 platform ([HG-U133_Plus_2] Affymetrix Human Genome U133
Plus 2.0 Array), [72]GSE1297 was built on [73]GPL96 Affymetrix Human
Genome U133A Array. All data were freely accessible online. The
[74]GSE5281 expression profile consists of 161 samples and
approximately 55,000 transcripts from 74 disease-free patients and 84
patients with AD. There are 31 samples in the [75]GSE1297 dataset,
including 9 normal and 22 AD subjects with different severities. The
correlation between each gene expression and clinical information of
the 31 subjects was verified by the Neurofibrillary Entanglement Score
(NFT) and the Minimental Examination (MMSE) ([76]Blalock et al., 2004;
[77]Liang et al., 2007). [78]GSE122063 is a gene expression profiling
in frontal and temporal cortices from 36 patients with vascular
dementia, 56 diseases (AD), and 44 non-demented controls (Control)
obtained from the University of Michigan. Brain Bank Mild cognitive
impairment (MCI) is a precursor to the development of AD. The
[79]GSE150693 data were downloaded from the GEO database from the
[80]GPL21263 platform, among them 197 miRNA samples from MCI sera
included 83 patients who were converted from MCIs to AD, and 114
patients who did not convert from MCI to AD ([81]Shigemizu et al.,
2020).
Application of weighted correlation network analysis algorithm
First, the [82]GSE1297 sample data obtained from the GEO database was
processed by using the log scale robust multi-array analysis (RMA).
Then the missing values in the samples were filled by using the
K-neighborhood algorithm to improve the accuracy and usability of the
data. Next, expression relationships among genes were measured by
correlation coefficients to construct weighted co-expression networks
in the absence of outlier samples through cluster analysis of gene
samples. Based on this, gene modules were identified by using a
clustering algorithm as well as a dynamic shear tree algorithm. The
calculated correlation coefficients between gene modules and clinical
performance traits were used as the association criteria. Using the
Pearson correlation test can analyze the relationship between module
eigengenes and clinical features, and select the modules for focused
mining. Among them, clinical information includes: Age, age at death;
NFT, neurofibrillary tangle count; Braak, Braak stage; MMSE, adjusted
Minimental Status Exam; PMI, postmortem interval. Finally, the key
modules were further mined to find the key core genes.
Functional classification and pathway enrichment of differentially expressed
genes
The above differentially expressed genes (DEGs) screened from the
weighted correlation network analysis (WGCNA) were submitted to gene
ontology (GO) function enrichment analysis, which was consisted of
cellular component (CC), biological process (BP), and molecular
function (MF), and was analyzed and visualized by using the R language
“cluster profile,” “GGploT2,” and “enrichPlot” packages. Kyoto
encyclopedia of genes and genomes (KEGG) pathway enrichment analysis of
DEGs in this study was performed by the Database for Annotation,
Visualization, and Integrated Discovery (DAVID) tools.^[83]2 KEGG
pathway visualization was analyzed by Sangerbox online pathway analysis
tool ([84]Shen et al., 2022). Enriched GO terms with adjusted value P <
0.01 and gene number> 12 were selected, and P > 0.05 of KEGG pathways
and gene numbers> 5 were considered statistically significant.
Analyzing differential expression
DEGs are processed using the “R” language “limma” package and calculate
adjusted P values and | logFC|. For [85]GSE5281 and [86]GSE97760 gene
expression profiles, Log2 Fold Change (FC) | > 1.0 and adjusted
P-values < 0.05 were selected as cut-off criteria.
Protein-protein interaction network and module analysis
WGCNA was used to directly set thresholds for gene saliency (GS) and
module feature (MM). Based on | MM| > 0.8, | GS| > 0.1 threshold value
for the standard to filter core genes, 269 key genes were selected from
789 genes in the “midnightblue” module. To explore the relationship
between these 269 hub genes and MAPT (AD is defined by the presence of
extracellular and intracellular neurofibrillary tangles composed of
hyperphosphorylated tau proteins) and AD disease ([87]Corbo and Alonso
Adel, 2011; [88]Bakota and Brandt, 2016; [89]Chong et al., 2018), the
Search Tool was used for the Retrieval of Interacting Genes (STRING
version 11.5) for known or predicted, direct (physical) and indirect
(functional) PPIs ([90]Szklarczyk et al., 2015). The interaction
confidence score of >0.4 was used as the cutoff value for statistical
significance. Subsequently, the PPI network was visualized by
Cytoscape3.9 software.^[91]3 To better search for proteins related to
Tau protein (MAPT), the algorithm between Centrality was adopted in
this study to find the important nodes [the calculation formula is (1)
and (2)].
We used CytoHubba, a plugin in Cytoscape, to compute nodes in the
protein network that efficiently convey data. In our study, the top 10
genes were identified as central genes.
[MATH: CB(Ni)=∑k,j<
/mi>=1gs
d(k,i,j),(k≠j)
:MATH]
(1)
sd(k, i, j) indicates that the shortest path from k to j passes through
i, which means, i is on the shortest path from k to j.
[MATH:
CB′(
Ni)=CB(Ni)∑k,j=1
g(k,j),(k≠j) :MATH]
(2)
The formula can be standardized to obtain (2), where the denominator
represents the number of paths between two points in the figure, that
is, the number of all paths.
Gene set enrichment analysis
For Gene Set Enrichment Analysis (GSEA), we obtained data from GSEA
([92]DOI: 10.1073/pans.0506580102)^[93]4 website for the GSEA software
(version 3.0), the samples were divided into high expression group
(≤50%) and low expression group (<50%) according to the expression
level of GRM5 and GRIA1. Obtained from the Molecular Signatures
Database ([94]DOI: 10.1093/bioinformatics/btr260),^[95]5 a subset of
“c3.mir.v7.4.symbols.gmt” and “c3.tft.v7.4.symbols.gmt” was downloaded
to evaluate related pathways and molecular mechanisms, based on gene
expression profile and phenotype grouping. The minimum gene set was 5.
The maximum gene set was 5,000, 1,000 resampling, P-value < 0.05 and
FDR < 0.25 were considered statistically significant, to verify the
binding relationship of miRNAs and transcription factors with the
target genes, and to indicate that they are on the relevant biological
pathway.
MicroRNAs associated with hub
Genes
The ENCORI (The Encyclopedia of RNA Interacts) database is used to
construct Hub genes and target miRNA interactions with them. We then
used the “R” software package “Survival” to integrate data on survival,
AD transformation status, and 35 features, and evaluated the prognostic
significance of these features in 197 samples by Cox and selection of
key microRNAs (miRNAs). Finally, these central genes and miRNA networks
were mapped using Cytoscape software.
The least absolute shrinkage and selection operator logistic regression
algorithm model was constructed
Here, we used the R package “glmnet” to integrate AD status and gene
expression data from [96]GSE5281 database for regression analysis using
least absolute shrinkage and selection operator (Lasso)-Cox. Tenfold
cross validation was performed to obtain the optimal model. We set the
Lambda value to 0.00208765767322578 and finally obtained four genes.
The model formula is: RiskScore = −0.00111259701091014 ×
TBP-0.00540566228324965 × CDK7 + 0.0013635456849925 × GRM5 −
0.000512469864405875 × GRIA1 and then, the database [97]GSE122063 was
used to verify the model to ensure the rationality of the model.
The key genes obtained from network analysis were verified with FpClass tool
The FpClass for predicting high confidence PPIs on a proteome-wide
scale, including proteins with few (low-degree proteins) or no known
binding partners (orphans). To identify co-expressed genes with the
queried genes, the tool considers two specific scores: the gene
expression score and the topological network score. The gene expression
score is based on the Pearson correlation of the gene expression
pattern. The topology score of the network shows whether genes exist in
the training data and the strength of their interaction ([98]Kotlyar et
al., 2015; [99]Sriroopreddy et al., 2019).
Drug selection of hub gene
Drugs targeting key genes were retrieved from the DrugBank database,
which is a web-enabled database.^[100]6 The database collects more than
7,800 drugs, including interactions and their targets, drug binding
data, drug-drug and drug-food interaction data. It also includes
hundreds of drug impact information on metabolite levels, protein
expression levels, hundreds of drug clinical trials and drug reuse
trials ([101]Wishart et al., 2018).
Results
Gene co-expression modules
To explore the co-expression patterns of mRNA in AD, we performed WGCNA
analysis on the [102]GSE1297 dataset. The dataset contained 22 columns
of AD samples and 9 columns of normal human samples, with a total of
22,283 genes. The raw data can be downloaded from the GEO database of
NCBI.^[103]7 The read-in data was first preprocessed by the rma
function in the Affy package in R language, which contains background
correction, normalization, calculation of expression to normalize the
data. Then, for the missing values in the data, we used the
K-neighborhood algorithm (KNN) to complement the missing values data.
Through this algorithm, genes with high similarity were identified to
supplement the null values. The absolute median difference (MAD) is the
median of the difference between each data in a set of data and the
median of that set of data. To filter out a large number of genes with
relatively constant expression in different gene samples, the top 7,200
genes of MAD values were selected for WGCNA analysis in this paper. To
explore the co-expression patterns of mRNA in AD, we performed WGCNA
analysis on the [104]GSE1297 dataset ([105]Figures 1D,F). To ensure
scale-free networks, we chose soft thresholds of b = 12 ([106]Figures
1A,B), used WGCNA packages as soft threshold power to generate
hierarchical clustering trees ([107]Figure 1C), and then we built
co-expression networks of associations between clinical features and
these modules. As shown in [108]Figure 2 the “Midnight” module of
[109]GSE1297 is significantly related to the clinical features of AD.
FIGURE 1.
[110]FIGURE 1
[111]Open in a new tab
Co-expression network of differentially expressed genes in AD. (A,B)
Soft threshold power selection. (C) Clustering tree dendrogram of
co-expression modules. Different colors represent distinct
co-expression modules. (D) Correlation analysis between midnightblue
module and clinical condition, each column represents one module,
individual rows represent clinical status. (E) Correlation between
modules. (F) Specific clinical information of the sample.
FIGURE 2.
[112]FIGURE 2
[113]Open in a new tab
Midnightblue module and clinical relevance. (A) Scatter plot between
module membership in midnightblue module and the clinical information
MMSE. (B) Scatter plot between module membership in midnightblue module
and clinical information ref number. (C,D) Scatter plot between module
membership in midnightblue module and clinical information NFT and
break, respectively.
Kyoto encyclopedia of genes and genomes pathway and gene ontology analysis of
differentially expressed genes
GO analysis of DEGs (including molecular function, biological processes
and cell composition analysis) and KEGG pathway analysis, respectively,
can provide valuable evidence of processes and pathways in which gene
sets may be involved. This type of information is crucial for
hypothesis development and for the design of future research.
Functional enrichment analyses of KEGG pathway enrichment and DEGs GO
function enrichment analysis for DEGs were performed using the DAVID
And R package “cluster profile” ([114]Tables 1, [115]2). The enriched
GO terms were divided into CC, BP, and MF ontologies. The results of GO
analysis indicated that DEGs were mainly enriched in BPs, including
modulation of chemical synaptic transmission, regulation of
trans-synaptic signaling, regulation of supramolecular fiber
organization, regulation of membrane potential, vesicle localization,
the establishment of organelle localization, protein-containing complex
disassembly and regulation of protein-containing complex assembly. MF
analysis showed that the DEGs were protein serine/threonine kinase
activity, tubulin binding, guanyl nucleotide binding, and geranyl
ribonucleotide binding. For the cell components, the DEGs were synaptic
membrane, presynapse, mitochondrial matrix, transporter complex,
neuron-to-neuron synapse, microtubule, and microtubule. In addition,
the results of KEGG pathway analysis showed that DEGs were mainly
enriched in pathways in retrograde endocannabinoid signaling, human
papillomavirus infection, tight junction, purine metabolism, and
dopaminergic synapse ([116]Figures 3, [117]4B and [118]Supplementary
Data 1).
TABLE 1.
Significantly enriched gene ontology (GO) terms of differentially
expressed genes (DEGs).
Category Term Description P-value Count
BP term GO:0051648 Vesicle localization 0.00000495 13
BP term GO:0051656 Establishment of organelle localization 0.0000306 17
BP term GO:0032984 Protein-containing complex disassembly 0.0000637 14
BP term GO:0043254 Regulation of protein-containing complex assembly
0.00015812 16
BP term GO:0050804 Modulation of chemical synaptic transmission
0.000179862 15
BP term GO:0099177 Regulation of trans-synaptic signaling 0.000184724
15
BP term GO:1902903 Regulation of supramolecular fiber organization
0.000237205 14
BP term GO:0042391 Regulation of membrane potential 0.00034946 15
CC term GO:0098793 Presynapse 0.0000018 19
CC term GO:0097060 Synaptic membrane 0.0000146 16
CC term GO:0005759 Mitochondrial matrix 0.0000279 18
CC term GO:1990351 Transporter complex 0.000311367 13
CC term GO:0098984 Neuron to neuron synapse 0.000348414 13
CC term GO:0005874 Microtubule 0.000822817 14
CC term GO:0043025 Neuronal cell body 0.002387617 14
MF term GO:0004674 Protein serine/threonine kinase activity 0.000143327
16
MF term GO:0015631 Tubulin binding 0.000968776 13
MF term GO:0019001 Guanyl nucleotide binding 0.001949535 13
MF term GO:0032561 Geranyl ribonucleotide binding 0.001949535 13
[119]Open in a new tab
TABLE 2.
Significant enrichment of the Kyoto encyclopedia of genes and genomes
(KEGG) pathway.
Category Term Description P-value Count
KEGG_PATHWAY hsa04723 Retrograde endocannabinoid signaling 0.004813647
8
KEGG_PATHWAY hsa05165 Human papillomavirus infection 0.018148818 11
KEGG_PATHWAY hsa04530 Tight junction 0.031975382 7
KEGG_PATHWAY hsa00230 Purine metabolism 0.034383145 6
KEGG_PATHWAY hsa04728 Dopaminergic synapse 0.038482344 6
[120]Open in a new tab
KEGG, Kyoto Encyclopedia of Genes and Genomes.
FIGURE 3.
[121]FIGURE 3
[122]Open in a new tab
Differentially expressed gene’s (DEG’s) GO analysis. (A) GO bioprocess
enrichment. (B) Cell composition enrichment results. (C) Results of
molecular function enrichment. The color of each bubble represents the
fitted p-value: the redder the color, the higher the concentration. (D)
The relationship between the classes. (E) The genes contained in the
more enriched category. GO, Gene Ontology.
FIGURE 4.
[123]FIGURE 4
[124]Open in a new tab
Kyoto encyclopedia of genes and genomes (KEGG) pathway map, TBP binding
site map, and key mRNA miRNAs prediction. (A) TBP transcription factor
binding site. (B) KEGG pathway diagrams of differentially expressed
genes (DEGs). (C) The miRNAs to GRM5 and GRIA1 bind are predicted. In
the network, thin lines represent sequence matches, green balls
represent miRNAs and blue balls represent mRNA. (D) The mRNA binding to
hsa-miR-425-5p and hsa-miR-186-5p was predicted. In the network, thin
lines represent sequence matching, green spheres represent mRNA, and
blue spheres represent miRNA. (E) The network diagram of Hub gene and
miRNA. The yellow ball represents the gene, the green represents the
miRNA, and the green triangle is the miRNA we finally screened.
Identification of differentially expressed genes in Alzheimer’s disease
We downloaded the series [125]GSE1297 and [126]GSE5281 datasets about
AD from the NCBI GEO database. After using the “limma” package in R
software, screening with the threshold of an adjusted p-value < 0.05
and | log2FC| > 1.0, 959 DEGs (320 upregulated and 639 downregulated)
were identified in the [127]GSE5281 dataset and DEGs were identified by
comparing AD samples with normal samples. 622 DEGs (408 upregulated and
215 downregulated) were identified in the [128]GSE1297 dataset, the
DEGs were identified by comparing “control” with the “Severe” status of
AD samples. The volcano plot and heatmap analyses were used to
visualize the DEGs of the two data sets shown in [129]Figures 5A,B,
[130]6A,B, respectively. In addition, a Venn diagram analysis was
performed to evaluate the common DEGs between [131]GSE5281, WGCNA-hub,
and [132]GSE1297. As presented in [133]Figure 6C, Intersecting, with 10
hub genes identified two overlapping DEGs (CDK7 and GRIA1).
FIGURE 5.
[134]FIGURE 5
[135]Open in a new tab
The volcano plot, heat map, and pathway’s bubble diagram of DEGs. (A) A
heat map of DEGs in [136]GSE1297. (B) Volcano plot of DEGs in
[137]GSE1297. (C) Pathway bubble diagram of DEGs screened by weighted
correlation network analysis (WGCNA). There are also two horizontal
dashed lines in the panel (B), representing log2FC at –1 and 1, The
vertical dashed line represents the adjusted p-value at 0.05, RAD21L1
and USP9Y represented the genes with the largest and smallest
difference, respectively, GRIA1 and CDK7 are the key genes in this
study. DEGs, differentially expressed genes.
FIGURE 6.
[138]FIGURE 6
[139]Open in a new tab
The volcano plot, heat map, and Venn diagram of DEGs. (A) A heat map of
DEGs in [140]GSE5281. (B) Volcano plot of DEGs in [141]GSE5281. (C)
Venn diagram of [142]GSE5281, [143]GSE1297, and WGCNA-hub. There are
also two horizontal dashed lines in the panel (B), representing log2FC
at –1 and 1, The vertical dashed line represents the adjusted p-value
at 0.05, ID3 and SST represented the genes with the largest and
smallest difference, respectively. GRIA1 and CDK7 are the key genes in
this study.
Building protein-protein interaction networks and screening 10 key genes
A total of 177 nodes and 343 edges were involved in the final PPI
network from 269 proteins using the STRING tool, as presented in
[144]Figure 7. The top 10 genes evaluated by the algorithm Between
Centrality were adopted in the final PPI network. The results showed
that the protein, a regulatory protein involved in mitosis (CCNB1), was
the most outstanding one, with a score = 10,430, followed by the
predominant excitatory neurotransmitter receptors in the mammalian
brain which are activated in a variety of normal neurophysiologic
processes (GRIA1; score = 4,243). The TFIID basal transcription factor
complex plays a major role in the initiation of RNA polymerase II (Pol
II)-dependent transcription. Other proteins are TATA-binding protein
(TBP; score = 4,035), the microtubule-associated protein, tau (MAPT;
score = 4,000), Protein Phosphatase 2 Regulatory Subunit Bdelta
(PPP2R2D; score = 3,479), Cyclin-Dependent Kinase 7, involved in cell
cycle control and RNA polymerase II-mediated RNA transcription (CDK7;
score = 3,158), USO1 Vesicle Transport Factor (USO1; score = 2,905), a
gene essential for cell survival and DNA repair (PRPF19; score =
2,809), The protein is a metabotropic glutamate receptor, whose
signaling activates a phosphatidylinositol-calcium second messenger
system. This protein may be involved in the regulation of neural
network activity and synaptic plasticity (GRM5; score = 2,800) and
Mitochondrial Ribosomal Protein S15 (MRPS15; score = 2,510). [145]Table
3 shows the top 10 hub genes.
FIGURE 7.
[146]FIGURE 7
[147]Open in a new tab
PPI network of hub genes, PPI network of genes constructed in the
STRING database. Then, visualization of a STRING-derived network of
molecular interactions in Cytoscape pathway visualization and analysis
software, the 10 genes in the middle are the key genes to be screened,
and the closer the color is to red, the higher the score.
TABLE 3.
The hub genes.
Rank Name Score
1 CCNB1 10429.51
2 GRIA1 4242.677
3 TBP 4035.28
4 MAPT 4000.436
5 PPP2R2D 3479.117
6 CDK7 3157.932
7 USO1 2904.839
8 PRPF19 2809.138
9 GRM5 2799.903
10 MRPS15 2510.792
[148]Open in a new tab
Prediction of microRNAs and identification of common differentially expressed
miRNAs
In addition, we explored the predicted miRNAs of GRM5 and CRIA1 in AD
patients using the ENCORI platform to establish potential miRNA
messenger RNA (mRNA) regulatory networks ([149]Wu et al., 2021). We
found 53 predicted differentially expressed microRNAs ([150]Figure 4C).
Also, 140 specifically dysregulated differentially expressed miRNAs
(DEmiRNAs) (58 downregulated and 82 upregulated DEmiRNAs) in blood
samples of patients with AD were identified and used in our work,
according to the report by [151]Leidinger et al. (2013;
[152]Supplementary Data 2). The DEmiRNAs from AD were intersected with
the 82 upregulated DEmiRNAs from blood samples; the common DEmiRNAs
were identified. Respectively, hsa-miR-186-5p and hsa-miR-425-5p were
identified to be potentially involved in the regulation of GRM5 and
GRIA1 ([153]Figures 8C,E).
FIGURE 8.
[154]FIGURE 8
[155]Open in a new tab
Selection of key miRNAs. (A) AD conversion curves of hsa-miR-186-5p.
(B) ROC analyses for the 365, 1,095, and 1,825 time points were
performed using the ROC function of the R package proc (version
1.17.0.1). (C) The intersection of predicted miRNAs and abnormal miRNAs
in blood. (D) hsa-miR-425-5p. (E) They are, respectively, the
complementary base sequence matching of hsa-miR-186-5p and target gene
GRM5, and the complementary base sequence matching of hsa-miR-425-5p
and target gene GRIA1.
Prediction of transcription factor binding to target genes
We downloaded the TBP’s primary structure of the amino acid sequence
from the Uniprot database,^[156]8 and then found the TBP transcription
factor binding site features in the JASPAR database^[157]9 ([158]Figure
4A). Then, 1,000 bp promoter regions of GRM5 and GRIA1 sequences were
downloaded from the UCSC database^[159]10 ([160]Supplementary Data 3).
Finally, their promoter sequences were entered into the JASPAR database
to predict regulatory transcription factors. Even when the threshold
was set as 80%, TBP could be found to be the transcription factor of
both. For a more accurate prediction of transcription factors, we used
the HumanTFDB database.^[161]11 Interestingly, not only TBP is a
transcription factor of GRM5 and GRIA1, but CDK7 is also a
transcription factor of both (detailed results are provided in the
[162]Supplementary Data 4).
Drug selection results
Fifteen drugs targeting the four core genes were selected based on drug
and target information contained in DrugBank database. Of these, nine
medications have been approved, six are investigational and
experimental. All four GRIA1 drugs are primarily responsible for the
activity of AMPA glutamate receptors and extracellular glutamate gated
ion channels. Isoflurane (DB00753), Methodflurane (DB01028), Desflurane
(DB01189), Sevoflurane (DB01236) are all known to be an antagonist
medication. Isoflurane is a general inhaled cosmetic used in surgery.
Methoxyflurane is a general inhalational anesthetic used for the
induction and maintenance of general anesthesia. Desflurane is a
general inhalational anesthetic for both inpatient and outpatient
surgery in adults. Sevoflurane is an inhalational anesthetic agent used
for the induction and maintenance of general anesthesia during surgical
interventions. Additional information on the features and specific
usage of other gene-targeted drugs is presented in [163]Supplementary
Data 6 and in [164]Table 4.
TABLE 4.
Fifteen drugs targeting key genes screened from DrugBank database.
genes Drug name and ID Status Type
TBP Quercetin (DB04216) Experimental, investigational
TBP Chloroquine (DB00608) Approved, investigational, vet_approved
Inhibitor
TBP Revusiran (DB16309) Investigational Regulator
CDK7 Phosphonothreonine (DB02482) Experimental
CDK7 Alvocidib (DB03496) Experimental, investigational
CDK7 SNS-032 (DB05969) Investigational
CDK7 Seliciclib (DB06195) Investigational
CDK7 Trilaciclib (DB15442) Approved, investigational Inhibitor
GRIA1 Isoflurane (DB00753) Approved, vet_approved Antagonist
GRIA1 Methoxyflurane (DB01028) Approved, investigational, vet_approved
Antagonist
GRIA1 Desflurane (DB01189) Approved Antagonist
GRIA1 Sevoflurane (DB01236) Approved, vet_approved Antagonist
GRM5 Imipramine (DB00458) Approved Inhibitor
GRM5 Disopyramide (DB00280) Approved Inhibitor
GRM5 Dalfampridine (DB06637) Approved Antagonist
[165]Open in a new tab
Discussion
Due to the heterogeneity of AD pathology, there is a lack of sufficient
efficacy in treating AD ([166]Lam et al., 2013; [167]Byun et al.,
2015). In the past decades, neurodegenerative therapy for AD has made
more and more progress, and drug resistance is often allowed in
traditional histology ([168]Nandigam, 2008). Hence, identifying more
appropriate molecular regulation of transcriptional clusters is
critical to guiding personalized AD treatment. In this work, the top
7,200 genes based on MAD values were chosen for WGCNA analysis. The
“Midnight” module of [169]GSE1297 was chosen as it is significantly
associated with clinical AD ([170]Figure 2). Subsequently, according to
the cytoHubba plug-in of Cytoscape, we made use of the algorithm
Between Centrality to screen 10 key DEGs as hub genes in the PPI
network, including CCNB1, GRIA1, TBP, MAPT, PPP2R2D, CDK7, USO1,
PRPF19, GRM5, and MRPS15. It’s suggesting that these genes may play
important role in the mechanism of AD and the specific flow chart is
shown in [171]Figure 9G.
FIGURE 9.
[172]FIGURE 9
[173]Open in a new tab
Protein-to-protein correlation graph ([174]GSE5281 data set). (A–C)
Respectively, the correlation map of GRIA1 and TBP, CDK7, MAPT. (D–F)
Respectively, the correlation map of GRM5 and TBP, CDK7, MAPT. (G) The
flow-process diagram.
Luckily, we discovered that in GO analysis, GRM5 and GRIA1 were
predominantly enhanced in the synaptic membrane of CC (cellular
component), Neuron to Neuron synapse, postsynaptic membrane,
postsynaptic density, and postsynaptic specialization. In BP
(biological process), GRM5 and GRIA1 were predominantly enriched in the
modulation of chemical synaptic transmission, regulation of
trans-synaptic signaling, and regulation of membrane potential. More
and more studies have shown that neuron synapse and AD are closely
related, synapse loss and Tau pathology are hallmarks of AD and other
tauopathies ([175]Yoshiyama et al., 2007; [176]Skaper et al., 2017;
[177]Dejanovic et al., 2018; [178]Kurucu et al., 2022). In KEGG pathway
analysis, GRM5 and GRIA1 are co-enriched in glutamatergic synapse and
retrograde endocannabinoid signaling pathways. TBP, GRM5, and GRIA1
were co-enriched in Huntington’s disease pathway. Both Huntington’s
disease (HD) and Alzheimer’s disease (AD) are similar, albeit
clinically distinct, neurodegenerative disorders that share pathologic
features related to selective brain injury ([179]Naia et al., 2017;
[180]Mukherjee, 2021). TBP and CDK7 are basal transcription factors
among them. In addition, through FpClass analysis, we know that TBP and
CDK7 have high reliability based on gene expression and topological
network analysis, and GRM5 and GRIA1 also have high scores
([181]Supplementary Data 7). In the verification in the GeneMANIA
database, they still have the same relationship network
([182]Supplementary Data 5, Figure A).
All three nuclear RNA polymerases have a subunit of the TATA binding
protein (TBP) complex required for transcription ([183]Davidson et al.,
2004). In all fields of life, transcription regulation by DNA-dependent
three-types RNA polymerase is largely achieved at the initial level.
Among them, TBP is the most conservative initiation factor, which is
crucial to the transcription initiation of all ancient eukaryotes, and
is the only initiation factor required for the complete start of RNA
polymerase in all eukaryotes ([184]Kramm et al., 2019; [185]Mishal and
Luna-Arias, 2022). Back in 2004, they demonstrated significant
differences in TBP between the AD group of normal subjects and
patients. In addition, their number and distribution are not
proportional to the number and the distribution of positive tau or
β-amyloid structures, and the TBP accumulates in AD brains, suggesting
that TBP might be a contributing factor to AD due to its own
entanglements ([186]Reid et al., 2004). Evidence implicating TBP in the
molecular mechanism of several neurodegenerative diseases has emerged
in the past few years. TBP may contribute to these diseases through a
loss of normal function (likely to be catastrophic to a cell) which
affects neurodegenerative diseases, such as AD and Huntington’s disease
([187]van Roon-Mom et al., 2005; [188]Bech et al., 2010).
The transcription cycle of RNA polymerase II is regulated by a group of
cyclin dependent kinases (CDK). CDK7 related to TFIIH, a transcription
initiation factor, is not only an effector of RNA polymerase II
phosphorylation and other targets in the transcription mechanism, but
also a CDK activated kinase involved in transcription, and also plays a
key role in regulating eukaryotic cell division ([189]Sansó and Fisher,
2013; [190]Fisher, 2019). There is growing evidence that CDK regulates
the transcriptional cycle of RNA polymerase II. Specific CDKs are
considered as important molecular mechanisms in the transcription
cycle, and describe that recently emerged transcriptional CDK can serve
as promising drug targets in cancer ([191]Fisher, 2017; [192]Parua and
Fisher, 2020). Moreover, some studies show that neuronal activity of
CDK7 in hippocampus is related to aging and AD ([193]Zhu et al., 2000).
CDK7 is also required for activity-Dependent expression of Neuronal
genes, synaptic plasticity at long-term, and memory at long-term
([194]He et al., 2017). It can be hypothesized that CDK7 is a promising
drug target for AD.
Although GRM5 (coding for metabotropic glutamate receptor 5, mGluR5) is
a promising target for treating cognitive deficits in schizophrenia and
AD, its association with cognitive and brain phenotypes within this
disorder has received little attention. Early researchers found that
the common genetic variation of GRM5 in schizophrenic patients would
affect cognitive function, hippocampal volume and hippocampal mGluR5
protein level compared with the healthy control group, among them, the
metabotropic glutamate receptor subtype 5 (mGluR5) is encoded by GRM5
gene, which is an attractive new drug target for the treatment of
schizophrenia ([195]Matosin et al., 2018). MGluR5 is a postsynaptic
G-protein coupled glutamate receptor, which is closely related to
several key cellular processes destroyed in schizophrenia. In
preclinical models of schizophrenia, positive mGluR5 modulators have
shown encouraging therapeutic potential, especially in the treatment of
cognitive dysfunctions ([196]Matosin et al., 2017). Early in the course
of AD inhibition, some have attempted to do so by downregulating
cholinergic receptors and glutamate receptors, which, as the disease
progresses, play key roles in inflammation and oxidative stress,
thereby participating in the neurodegenerative process. Accumulating
evidence suggests that perturbation of the excitatory amino acid
L-glutamate (L-Glu) in systems may underlie the pathogenesis of hypoxia
ischemia, epilepsy, and chronic neurodegenerative diseases (e.g.,
Huntington’s disease and AD) ([197]Hynd et al., 2004; [198]Schaeffer
and Gattaz, 2008). It is reported that glutamate excitatory
neurotransmission is an important process in learning and memory, which
is seriously damaged in AD, possibly due to β Amyloid peptides (1–42)
increase in relation to oxidative stress. The researchers also found
that some glutamate receptors are overactivated in AD. This sustained
mild activation may lead to neuronal damage and impaired synaptic
plasticity (learning). More and more evidence shows that glutamate
mediated neurotoxicity is involved in the pathogenesis of AD, and this
metabolic glutamate receptor (GRM5), whose signal transduction will
activate the second messenger system of calcium phosphatidyl inositol,
will also participate in the regulation of neural network activity and
synaptic plasticity, which makes us speculate that GRM5 may be a very
promising target for exploring treatment and improving AD disease
([199]Danysz and Parsons, 2003; [200]Doraiswamy, 2003; [201]Tanović and
Alfaro, 2006; [202]Świetlik et al., 2022).
GluA1 (also known as GluRA or GluR1) is a glutamate receptor subunit of
AMPA encoded by the GRIA1 gene, there is genome-wide association
between GRIA1 and schizophrenia ([203]Barkus et al., 2014; [204]Ang et
al., 2018). Hippocampal proteomics study has pointed out that proteins
involved in neuronal excitability and synaptic plasticity (e.g., GRIA1,
GRM3, and SYN1) were altered in both “normal” aging and AD ([205]Neuner
et al., 2017). Similarly, a bioinformatics study identified 464
differentially expressed genes that were modulated by silent
transcription factor (REST) between AD patients and controls. REST is
strongly associated with glutamatergic synapses and long-term
potentiation, and among them, GRIA1 shows a significant difference in
its tendency to change with REST ([206]Xu et al., 2021). From this, we
speculate GRIA1 also may be regulated by other transcription factors,
thereby affecting the progression of AD disease. From RNA-Seq
Expression Data from GTEx (53 Tissues, 570 Donors), the highest median
expression of GRM5 and GRIA1 all in Brain—Cerebellum
([207]Supplementary Data 5, Figure B).
A microRNA (miRNA) is a short non-coding RNA with regulatory functions
in a variety of biological processes, which has been implicated in many
cellular processes including cell proliferation, apoptosis, gene
expression, cellular differentiation, and development ([208]Krol et
al., 2010). Many studies have shown that miRNA often binds to target
mRNA through complementary sequences and induces translational
inhibition or target degradation as a negative regulation of mRNA
expression ([209]Towler et al., 2015). We tried to use the ENCORI
platform to predict the miRNAs that regulate GRM5 and GRIA1,
respectively, and then intersected with 82 up-regulated DEmiRNAs in the
blood samples of the previous experiment to determine the hub miRNAs
with significant differences. Respectively, hsa-miR-186-5p and
hsa-miR-425-5p were identified to be potentially involved in the
regulation of GRM5 and GRIA1. Among them, As shown in [210]Figures
8A,B,D, from the AD transformation state diagram and ROC curve, we
identified two significantly up-regulated miRNAs, hsa-miR-186-5p and
hsa-miR-425-5p may with good prognostic values. Recent studies have
found that hsa-miR-425-5p and hsa-miR-186-5p can be used as therapeutic
targets for other diseases. For example, hsa-miR-425-5p may promote
tumor occurrence and metastasis by activating CTNND1 related pathway
([211]Liu et al., 2020), hsa-miR-186-5p regulates TGFβ signaling
pathway through expression suppression of SMAD7 and SMAD6 genes in
colorectal cancer ([212]Bayat et al., 2021). In our study, GRM5 and
GRIA1 were found to be significantly downregulated at mRNA level and
regulated by the up-regulated hsa-miR-186-5p and hsa-miR-425-5p,
respectively ([213]Figures 4D,E). It has been reported that negatively
regulated miRNA-mRNA pairs contribute to the improvement of AD and
provide new reliable targets for the treatment of AD ([214]Fu et al.,
2019; [215]Qian et al., 2019; [216]Chen et al., 2021).
Anomalous aggregation of microtubule associated tau protein (MAPT) is a
prominent pathological feature in various neurodegenerative diseases
including AD ([217]Zhou and Wang, 2017; [218]Strang et al., 2019). From
in [219]Figure 9, we can observe that the correlation coefficients of
GRM5 and GRIA1 with MAPT were 0.69 and 0.18, respectively, with P <
0.05. A reliable index of AD-related cognitive state at a particular
time point is the MiniMental State Examination (MMSE) ([220]Clark et
al., 1999). Furthermore, dementia-NFT of the senile type has been
termed tangle-only dementia, NFT-predominant form of SD and limbic NFT
dementia ([221]Hyman, 1997). Therefore, MMSE and NFT were chosen as our
main markers to quantify AD progression ([222]Schmitt et al., 2000;
[223]Yamada, 2003). Since NFT scores increase and MMSE scores decrease
with AD severity, genes upregulated with AD could only correlate
positively with NFT or negatively with MMSE, whereas genes
downregulated with AD could only positively correlate with MMES or
negatively correlated with NTF scores. In our study ([224]Figure 10),
the significantly down-regulated gene GRM5 was positively correlated
with MMSE with a correlation coefficient of 0.58 and negatively
correlated with NFT with a correlation coefficient of −0.4, GRIA1 and
MMSE were also positively correlated with a correlation coefficient of
0.52 and negatively correlated with NFT with a correlation coefficient
of −0.5. At the same time, their p < 0.05. As shown in [225]Figures
9A–F, we also can find that the transcription factors TBP and CDK7 are
well correlated with their downstream regulated genes. For the
transcription factors TBP and CDK7, and the differential expression
changes of their regulated target genes GRM5 and GRIA1. It can be seen
from [226]Figures 5, [227]6 that the expression of CDK7 and GRIA1 in AD
group is two times lower than that in normal group in two databases.
The expression of TBP and GRM5 also decreased significantly compared
with the normal group. Grouped according to gene expression profiles
and phenotypes, related pathways and molecular mechanisms were
evaluated by using GSEA software to validate the regulate of
hsa-miR-425-5p, hsa-miR-186, and TBP to target genes GRIA1 and GRM5,
respectively ([228]Figure 11).
FIGURE 10.
[229]FIGURE 10
[230]Open in a new tab
Boxplots of hub mRNA and correlation with clinical information. (A–D)
Based on the difference boxplot of TBP, GRIA1, GRM5, and CDK7 between
the normal group and AD group ([231]GSE5281 data set). (E–H) The
correlation between GRIA1, GRM5, and clinical information MMSE, NFT,
respectively.
FIGURE 11.
[232]FIGURE 11
[233]Open in a new tab
GSEA software was used to assess miRNA and TBP in similar biological
pathways with target genes. (A,B) Grouped according to gene expression
profiles and phenotypes, and evaluated related pathways and molecular
mechanisms to validate the regulate of the hsa-miR-425-5p and TBP to
the target gene GRIA1. (C,D) Validated regulate of hsa-miR-186 and TBP
to the target gene GRM5, and indicate that they are on the relevant
biological pathway. (E,F) They are GRIA1 and GRM5, respectively, and
their corresponding miRNA-related biological pathways.
We also performed regression analysis using Lasso-Cox to construct and
validate a prognostic signature associated with the four key genes,
which had significant predictive value for AD patients and achieved
high AUC values close to 90% in both databases ([234]Figure 12). In
addition, we screened 15 drugs targeting key genes from the DrugBank
database, most of which are used to relieve pain, systemic anesthetics,
inhibit signaling, inhibit cyclin-dependent kinases, inhibit multiple
enzyme targets (including CDK7) and change the growth phase of treated
cells, treat depression, it is used in potassium channel blockers,
ameliorates multiple sclerosis (MS), and so on. We can try to study
these drugs, and maybe they may help us to alleviate and treat AD
([235]Supplementary Data 7). In conclusion, based on the bioinformatics
analysis in this study, we propose GRM5 and GRIA1 as novel potential
prognostic markers of AD and suggest miRNAs and transcription factors
that may be related to AD and regulate GRM5 and GRIA1. We hope our
findings will inform future research to improve outcomes for patients
and try to mitigate and treat AD.
FIGURE 12.
[236]FIGURE 12
[237]Open in a new tab
Least absolute shrinkage and selection operator (LASSO) regression
algorithm and model prognostic evaluation. (A,B) The model based on the
[238]GSE5281 dataset with four AD-related genes was determined by
Lasso-Cox algorithm. (C) Heatmap of risk score distribution and
prognostic 4 gene signatures in the [239]GSE5281 database. (D,E) ROC
curve and corresponding AUC value of two databases. (D) Database for
model building [240]GSE5281. (E) Database for validation
[241]GSE122063.
Conclusion
Our research identified a miRNA/TF-gene network that is potentially
relevant for AD. The four hub genes (including TBP, CDK7 GRIA1, and
GRM5) were markedly downregulated, which may have a critical influence
on the pathophysiological mechanism of AD. Two potential target miRNAs
(hsa-miR-425-5p and hsa-miR-186-5p) were furthermore forecast. The
significantly down-regulated transcription factors TBP and CDK7 may be
involved in the transcription of GRM5 and GRIA1, thereby affecting
their expression and thus the progression of AD. TBP accumulating in
the AD brain, localizing to neurofibrillary tangle structures, may be a
contributing factor in AD. CDK7 is necessary for Activity-Dependent
Neuronal Gene Expression, Long-Term Synaptic Plasticity, and Long-Term
Memory, some studies have demonstrated that CDK7 neuronal activity in
the hippocampus is linked to aging and AD. The downregulation of both
TBP and CDK7 makes GRM5 and GRIA1 possible to cause neuronal damage and
impaired synaptic plasticity in the pathogenesis of AD. As the
metabotropic glutamate receptor encoded by GRM5 gene and the AMPA
glutamate receptor subunit encoded by GRIA1 participate in the
perturbed glutamatergic neurotransmission. Meanwhile, up-regulated
hsa-miR-425-5p and hsa-miR-186-5p inhibited the transcription of GRIA1
and GRM5 by binding to them. These findings may contribute to the early
diagnostic strategies, treatment targets and prognostic markers of AD
disease.
Data availability statement
The original contributions presented in this study are included in the
article/[242]Supplementary material, further inquiries can be directed
to the corresponding authors.
Author contributions
QZ and PY designed the research study and wrote the original article.
QZ, WG, and PY performed the data analysis. QZ, XP, YS, YW, and CP
participated in the conception of the study and revision of the
manuscript. YW was responsible for guiding the research framework,
including concepts and theories. All authors approved the final version
to be published.
Footnotes
^1
[243]https://www.ncbi.nlm.nih.gov/geo/browse/
^2
[244]https://david.ncifcrf.gov/
^3
[245]www.cytoscape.org/
^4
[246]http://software.broadinstitute.org/gsea/index.jsp
^5
[247]http://www.gsea-msigdb.org/gsea/downloads.jsp
^6
[248]www.DrugBank.ca
^7
[249]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE1297
^8
[250]https://www.uniprot.org/
^9
[251]https://jaspar.genereg.net/
^10
[252]http://genome.ucsc.edu/
^11
[253]http://bioinfo.life.hust.edu.cn/HumanTFDB#!/
Funding
This work was funded by the Sichuan Normal University. This work was
supported in part by the National Natural Science Foundation of China
(Grant No. 61988102).
Conflict of interest
The authors declare that the research was conducted in the absence of
any commercial or financial relationships that could be construed as a
potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors
and do not necessarily represent those of their affiliated
organizations, or those of the publisher, the editors and the
reviewers. Any product that may be evaluated in this article, or claim
that may be made by its manufacturer, is not guaranteed or endorsed by
the publisher.
Supplementary material
The Supplementary Material for this article can be found online at:
[254]https://www.frontiersin.org/articles/10.3389/fnagi.2022.1069606/fu
ll#supplementary-material
[255]Click here for additional data file.^ (8.6MB, zip)
References