Construction of a ceRNA immunoregulatory network related to the development of vascular dementia through a weighted gene coexpression network analysis.

BACKGROUND: To date, vascular dementia (VaD) lacks effective treatment in clinical practice. There is also growing evidence that VaD may be closely related to the immune response. The development of high-throughput technology, and the recently discovered group of new mediators called competitive endogenous RNAs (ceRNA), provides a unique opportunity to study the immunomodulation of VaD.
METHODS: In this study, we used gene expression profiles in the Gene Expression Omnibus (GEO) database to obtain immune-related gene coexpression modules through a weighted gene coexpression network analysis (WGCNA) and gene enrichment analysis. We extracted and merged long non-coding RNA (lncRNA) and microRNA (miRNA) expressions from the GEO database and mapped them with related databases. Subsequently, we used Cytoscape to construct a lncRNA-miRNA-mRNA network, and then we performed an enrichment analysis on the mRNAs in the network to determine their regulatory function. Subsequently, we used an ImmuCellAI immune infiltration analysis and constructed a ceRNA sub-network of related immune cells. Finally, we conducted a gene set enrichment analysis (GSEA) to determine the potential regulatory pathways of the key factors.
RESULTS: As a result, we identified the blue module as the key module of immunity and constructed the related CeRNA network. Immune infiltration analysis showed that natural killer T (NKT) cells may be the key immune cells of VaD. Using a Pearson correlation analysis, we identified that B4GALT1, PPP1R3B, MICB, HHAT, DSC2, DNA2, SCARA3, and lncRNA NEAT1 may be the key factors of VaD. Our subsequent GSEA analysis showed that lncRNA NEAT1 may be regulated by NK cells and toll-like receptors.
CONCLUSIONS: Our research provides new therapeutic targets for vascular dementia from the immunological perspective for the first time, including B4GALT1, PPP1R3B, MICB, HHAT, DSC2, DNA2, SCARA3, and lncRNA NEAT1, and our research hopes to provide new treatment options for VaD. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Vascular dementia (VaD) is the second most common type of dementia. It is related to vascular disease, which leads to dysfunction and damage to the cerebrovascular system, including interruption of cerebral blood flow (CBF), brain damage, and eventual cognitive impairment and memory loss (1). In North America and Europe (2), VaD accounts for approximately 15% to 20% of dementia cases, while in Asia and developing countries this proportion has risen to approximately 30% (3). Essentially, VaD is caused by a decrease in the blood supply to the brain (4) and is due to partial or complete occlusion of cerebral perfusion, which is crucial for memory and cognitive processes. VaD reaches its peak in ischemic brain injury (5). In addition to vascular problems, risk factors such as diabetes, age, smoking, and alcohol cannot be ignored (6). Improving vascular endothelial dysfunction and restoring the increase in oxidative stress and neuroinflammation caused by hypoperfusion may provide a potential treatment strategy for preventing the cognitive dysfunction of VaD. From the morphological and structural point of view, the brain of VaD is mainly manifested as a tendency of brain atrophy, and VaD can lead to the decline of brain function, especially the comprehensive decline of intelligence-related functions, and the syndrome will decline to a certain degree. The current treatment for VaD is symptomatic rather than radical, as the number of known drug targets is insufficient to develop VaD therapies. Currently, no drug has been approved by the US FAD for VaD.

In the different forms of central nervous system (CNS) diseases, neuroimmune disorders are common phenomena (7). Recent research indicates that the immune system plays a vital role in VaD, and immune dysfunction has been shown to exist in VaD (8). The CD4+ CD25+ FoxP3 regulatory T cells present in VaD patients are significantly reduced (9), and monocytes and natural killer (NK) cells are reduced in VaD, but not in Alzheimer‘s disease (AD) (10). In addition, the immune response has recently been shown to be an important factor in promoting the development of cerebral small vessel disease, which also affects the progression of VaD (11).

High-throughput technology has recently gained attention, providing an excellent opportunity to identify VaD immune-related biomarkers (12). Long non-coding RNA (lncRNA) is the largest category of ncRNA with a length >200 nucleotides, but it is no significant protein coding potential (13). Although they are not involved in protein coding, growing evidence shows that lncRNAs have many biological and molecular functions, such as epigenetic regulation, signal transduction, cell differentiation, and immune regulation (14). MicroRNAs (miRNAs, a small family of non-coding RNAs that are important post-transcriptional regulators of gene expression) and lncRNAs are gene regulatory markers that regulate the translation of various proteins and biological processes in health and disease states and play an important role in neurological diseases (15). In 2011, the competitive endogenous RNA (ceRNA) hypothesis proposed attracted much attention (16). In this hypothesis, lncRNA transcripts act as ceRNAs or natural miRNA sponges. They communicate with each other and coregulate mRNAs by competing with shared miRNAs (17). Surprisingly, the amount of research on the association between lncRNAs and VaD is sparse (15).

In this study, we revealed for the first time the molecular mechanism of ceRNA immunoregulation during the progression of VaD. We constructed a VaD-related ceRNA network through a weighted gene coexpression analysis (WGCNA) method to ensure that the nodes were highly interconnected and that the network was highly reliable (18). Then, we performed an immune cell infiltration analysis, constructed a related ceRNA subnet and screened the key factors. Our findings may assist in screening VaD biomarker candidates and therapeutic targets. Our research process is shown in .

Figure 1

Flow chart of this research.

Methods

Microarray data

From the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo), which is a public functional genomics data repository, we retrieved and downloaded the VaD microarray expression profile data set (GSE122063 and GSE120584). GSE122063 data, which is the lncRNA and mRNA data set, was collected by researchers at the University of Michigan from the frontal and temporal cortex of VaD, AD, and non-dementia control groups. The total sample size was 136 cases, including 36 cases of VaD and 44 control cases. The data of GSE120584, which is the miRNA data set, was collected from serum samples by researchers at the National Center for Geriatrics in Japan. The sample size was 1,601 cases, including 91 cases of AD, 56 cases of VaD, and 288 control cases. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Data preprocessing and the differential analysis

We downloaded GPL16699 and GPL21263 files from GEO. We transformed the probe into the corresponding gene symbol according to the annotation information. We used the “GEO2R” online analysis tool in the GEO database to analyze the difference between the two data sets and normalized them. For the mRNA in the chip, we set the threshold as an adjusted P value <0.05, and for the lncRNA and miRNA in the chip, the threshold was |log2FC| >0.5 and adjusted P value <0.05 (19).

The weighted gene coexpression network

We imported the gene expression matrix of the relevant VaD and control groups into the Image GP tool (http://www.ehbio.com/Cloud_Platform/front/#/) and constructed a coexpression network of the VaD and control groups’ differentially expressed genes through a WGCNA algorithm. We utilized WGCNA to construct a coexpression module, and screen disease candidate targets. The Pearson correlation coefficient was used to calculate the linear correlation between each gene pair. We utilized the power function to transform the matrix into a signed adjacency matrix, built a scale-free network, and then transformed it into a topological overlap matrix (TOM). Then, we performed hierarchical clustering of the highly co-expressed genes. Next, we used the dynamic tree cutting algorithm to cut the clustering tree branches to generate modules. The first principal component is the module eigengene (ME), which summarizes the gene expression overall level of each module. MEs are used to assess module-trait associations. Modular Significance (MS) and Gene Significance (GS) are used to calculate the expression patterns of modules related to disease traits. In this study, the dynamic tree cutting algorithm can be used to detect the module, and the minimum module size was 30 (20).

Module enrichment analysis

We used the Enrichr database (https://maayanlab.cloud/Enrichr/) to analyze each module of the WGCNA results using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis to study the function of related mRNAs. We then screened out the immune-related modules for the ceRNA network construction.

Construction of the lncRNA-miRNA-mRNA ceRNA network

According to the theory that lncRNAs act as miRNA sponges to influence miRNAs and further regulate mRNA expression, we constructed a ceRNA network. We used the starBase database (http://starbase.sysu.edu.cn/index.php) and the miRTarbase database (http://mirtarbase.mbc.nctu.edu.tw/php/index.php) (21) to search for mRNAs targeted by differential miRNAs. In order to improve the validity of the search results, we included miRNA-targeted mRNAs in all two databases, and then we cross-referenced these mRNAs with the mRNAs of the immune module and performed a correlation test to confirm that the miRNA and mRNA expressions were negatively correlated. We used the starBase (http://starbase.sysu.edu.cn/index.php) database to predict the lncRNAs regulated by differential miRNAs by taking the intersection with the differential lncRNAs obtained through the gene chip analysis, and then we performed a correlation test to confirm that the miRNA and mRNA expressions were negatively correlated. In the above process, we used the Cytoscape 3.7.2 visualization to construct a ceRNA network diagram.

mRNAs of ceRNA network Enrichment analysis

We used Enrichr to perform the GO and KEGG pathway enrichment analysis on mRNAs in this ceRNA network.

Immune infiltration analysis

ImmuCellAI (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/) is a method based on gene set characteristics and is used to accurately appraise the abundance of 18 T cell subgroups and 6 other important immune cells from gene microarray, including B cells, NK cells, and monocyte, macrophage, neutrophil and dendritic cells (DCs) (22). We used this tool to investigate the difference in immune cell infiltration between the VaD and control expression groups.

Construction of the macrophage-related ceRNA subnetwork

Pearson’s method was used to analyze the correlation coefficient between the value of the NKT cell infiltration and the expression value of mRNAs and lncRNAs (|R|>0.4, P<0.05) in the ceRNA network. Then we analyzed the correlation between the expression values of the lncRNAs and mRNAs in the results, and finally screened the NKT cell ceRNA sub-network construction with an R value >0.4.

Determination of the regulatory pathways of the key factors

After grouping according to the median expression value of key lncRNAs, we utilized gene set enrichment analysis (GSEA) to analysis the KEGG pathway and GO process of the gene microarray. The standardized enrichment score (NES) was acquired and utilized to identify the key pathways.

Statistical analysis

In this study, we used the online analysis tool GEO2R of the GEO database to analyze the difference of related chips, and used WGCNA to analyze the mRNAs matrix. Subsequently, we applied enrichment analysis to related mRNAs, and performed ImmuCellAI immune infiltration analysis on mRNAs. Finally, we used Pearson correlation analysis to analyze the ceRNA network, and finally performed GSEA analysis on the mRNA matrix.

Results

Differentially expressed mRNAs, lncRNAs, and miRNAs in VaD

We used the online analysis tool GEO2R of the GEO database to analyze the difference between the GSE122063 and GSE120584 data sets. After batch correction, we obtained 6,876 VaD-related differential mRNAs using an adjusted P value <0.05, of which 3,178 were upregulated and 3,698 were downregulated. In addition, we also obtained 94 differential miRNAs based on an adjusted P value <0.05 and |log2FC| >0.5, of which 89 were upregulated and 5 were downregulated (the related heat map and volcano map are shown in ). We also obtained 1,895 differential lncRNAs, of which 925 were upregulated and 970 were downregulated (the related heat map and volcano map are shown in ).

Figure 2

Volcano map and cluster map. (A) Volcano map of miRNAs; (B) cluster map of miRNAs; (C) volcano map of lncRNAs; (D) cluster map of lncRNAs.

Coexpression module construction and immune module identification

We used the WGCNA algorithm to construct a coexpression network and modules for the VaD and control group samples. Based on the scale-free topological structure of R2=0.85, the Pearson correlation matrix of the gene was transformed into a strengthened adjacency matrix according to the power r=6 (). All selected genes were clustered using the dissimilarity measurement method based on the topological overlap matrix (TOM), and the tree was divided into six modules with different colors based on the dynamic tree cutting algorithm (). Then we calculated and plotted the relationship between each module and its corresponding clinical features (). We used the Enrichr database to perform the enrichment analysis of KEGG and GO for each module. The KEGG pathway mainly included staphylococcus aureus infection; leishmaniasis; inflammatory bowel disease (IBD); phagosome; and tuberculosis (). It can be seen that the blue module may be the key module of the immune response, and the blue module contains a total of 393 mRNAs. Through enrichment analysis, we found that the GO biological process of the blue module mainly included the cytokine-mediated signaling pathway (GO:0019221); neutrophil activation involved in immune response (GO:0002283); neutrophil degranulation (GO:0043312); neutrophil-mediated immunity (GO:0002446); and inflammatory response (GO:0006954) ().

Figure 3

Weighted co-expression network and gene module analysis. (A and B) Soft threshold selection process; (C) clustering dendrogram: each color represents a specific co-expression module. In the colored rows below the tree diagram, the two colored rows represent the original module and the merged module respectively. (D) Heat map of the correlation between vascular dementia (VaD) and the characteristic genes of the module.

Figure 4

Module enrichment analysis and Venn diagram. (A) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis result of blue module; (B) Gene Ontology (GO) analysis result of blue module; (C) Venn plot of mRNAs predicted by differentially expressed miRNAs and blue modules; (D) Venn plot of lncRNAs predicted by differentially expressed miRNAs.

We predicted the mRNAs regulated by differential miRNAs through the starBase and miRTarbase databases, including 8,001 in the starBase database and 6,336 in the miRTarbase database. We selected the mRNAs contained in the two databases and intersected the immune-related modules, and a total of 75 mRNAs were obtained (). In addition, we obtained 46 mRNAs after correlation analysis. We predicted the miRNA-related lncRNAs through the starBase database, and obtained a total of 983 lncRNAs, and then intersected these with the differential lncRNAs to get a total of 16 lncRNAs (). A total of nine lncRNAs were obtained after correlation analysis. We used Cytoscape 3.7.2 to construct the relevant ceRNA network ().

Figure 5

Enrichment analysis of the immune-related ceRNA networks and mRNAs. (A) ceRNA network; (B) Gene Ontology (GO) enrichment analysis of mRNAs in the ceRNA network; (C) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of mRNAs in the ceRNA network.

Enrichment analysis

We used the Enrichr database to perform the enrichment analysis of GO and KEGG for the mRNAs. Through enrichment analysis, we found that the GO biological process of the blue module mainly included regulation of smooth muscle cell proliferation (GO:0048660); positive regulation of cell proliferation (GO:0008284); positive regulation of smooth muscle cell proliferation (GO:0048661); regulation of establishment of endothelial barrier (GO:1903140); and positive regulation of granulocyte macrophage colony-stimulating factor production (GO:0032725) (). The KEGG pathway mainly included the NF-kappa B signaling pathway; Epstein-Barr virus infection; TNF signaling pathway; legionellosis; and the cytosolic DNA-sensing pathway ().

Immune cell abundance analysis

ImmuCellAI was used to estimate the abundance of 24 immune cells in the GSE cohort based on microarray data. Among them, CD4_naive, natural regulatory T (nTreg), iTreg, Th1, natural killer T (NKT) cells, mucosal-associated invariant T cells (MAIT), B cells, monocytes, macrophages, natural killer (NK) cells, neutrophils, and infiltration scores were significantly different between the VaD and control groups. Of these, the difference in NKT cells was the most significant, with the VaD group showing significantly less abundance than the control group (P=2.2e-8) ().

Figure 6

Immune infiltration analysis and the correlation between lncRNAs and mRNAs. (A) ImmuCellAI immune infiltration analysis; (B) correlation analysis of lncRNA and mRNA in the ceRNA sub-network of natural killer T (NKT) cells.

Construction of the NKT cell-related ceRNA sub-network

We used GSE122063 to identify key genes related to the ImmuCellAI infiltration score of the NKT cells. After correlation analysis, it was found that one lncRNA and 14 mRNA had a strong correlation with NKT cell infiltration (https://cdn.amegroups.cn/static/public/atm-21-1717-1.zip). The subsequent correlation analysis established the NKT cell ceRNA sub-network constructed with B4GALT1, PPP1R3B, MICB, HHAT, DSC2, DNA2, SCARA3, and lncRNA NEAT1 (), and indicated that lncRNA NEAT1 may be the key lncRNA identified in our research.

GSEA analysis

We applied GSEA analysis to GSE122063 to clarify the regulatory pathways of the key prognostic factors through median grouping of lncRNA NEAT1 (). It can be seen that these key factors are closely related to immunity, and surprisingly, compared with the lncRNA low expression group, we found that the significantly enriched signaling pathways included the NATURAL_KILLER_CELL, and TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY (), and GO analysis also shows that lncRNA is closely related to the immune process ().

Figure 7

The ceRNA network is related to the natural killer T (NKT) cell infiltration score. (A) The NKT cell-related ceRNA sub-network; (B) the Gene Set Enrichment Analysis (GSEA) KEGG in response to lncRNA NEAT1; (C) the GSEA GO in response to lncRNA NEAT1.

Discussion

VaD refers to a syndrome associated with clinical stroke or subclinical cerebrovascular injury, which affects at least one cognitive function. The disease is a heterogeneous encephalopathy, accounting for more than 20% of dementia cases (23). Studies have shown that effective treatment of VaD is more challenging than AD. The number of known drug targets is insufficient to develop VaD therapies (4). Therefore, the identification of new therapeutic factors is especially meaningful.

Immune regulation plays an important role in VaD. There is no research on VaD related to NKT cells. However, previous studies have shown that the mechanisms mediated by NKT cells are pro-inflammatory and immune regulatory, including the protective effect of B cell regulation. The production of sex antibodies, and the destructive effects of suppressing autoimmunity through cytokine production, DC crosstalk, and T-/B-cell interaction (24) is consistent with our immune infiltration results. In this study, we found a CeRNA immunoregulatory network in VaD. Then we established a CeRNA subnet related to NKT cells and identified several key factors, including B4GALT1, PPP1R3B, MICB, HHAT, DSC2, DNA2, and lncRNA NEAT1, which could be used as potential therapeutic targets for VaD.

As far as we know, this is the first study using a WGCNA method to construct a lncRNA-related ceRNA network in VaD, to construct its immune-related ceRNA, and to screen the key factors. We first determined the immune-related modules in the WGCNA results through enrichment analysis, and constructed the immune response-related ceRNA network, and then enriched the mRNAs in the network to confirm the regulatory function of the network. In the established NKT cell-related ceRNA sub-network, we found that lncRNA NEAT1 may be the key lncRNA for immune regulation. The research on human NEAT1 mainly focuses on its oncogene role in various cancers (25). One study has shown that lncRNA NEAT1 is an epigenetic inhibitor of hippocampal-dependent long-term memory formation and suggests a potential relationship with neuroinflammation (26). Another study showed that the expression of lncRNA NEAT1 was significantly upregulated in patients with ischemic stroke and aggravated cerebral ischemia-reperfusion (27). In addition, knockdown of lncRNA NEAT1 inhibits the activity of the pro-inflammatory AKT/STAT3 pathway, and significantly reduces the differentiation and apoptosis of microglia to the pro-inflammatory (M1) phenotype (28). Furthermore, the sponge mechanism of miR-24-3p can exacerbate cerebral ischemia-reperfusion injury (29).

Most of the mRNAs in the ceRNA subnetwork associated with NKT cells may also play an important role in VaD. B4GALT1 is a member of the b-1 4-galactosyltransferase gene (B4GALT) family and plays an important role in many biological events, including morphogenesis, mammalian fertilization, brain development, and cell adhesion (30). This gene may also play an important role in thrombosis. A recent study reported that the expression and activity of B4GALT1 mRNA in platelets were significantly increased, indicating that the gene was upregulated during late thrombosis and platelet release, thus positioning the gene as a regulatory gene for late thrombosis, possibly at the platelet release level (31). This suggests that the formation of thrombus may have an important influence on the progression of VaD. The protein phosphatase 1 regulatory subunit 3B gene (PPP1R3B) is located on chromosome 8p. The previous two GWASs indicate that PPP1R3B polymorphism is related to blood lipid traits. In animal studies, a high expression of PPP1R3B can affect blood lipid traits. In addition, it was found that PPP1R3B variants affect cardiovascular events and atherosclerosis (32), and these factors may cause vascular disease, which in turn may cause VaD. There are relatively few studies on the MICB gene related to VaD. The human genome encodes eight functional “stress-inducing ligands”: MICA, MICB, and ULBP1-6. All of these are recognized by a single receptor NKG2D, which is expressed on natural killer (NK) cells, cytotoxic T cells, and other T cell subpopulations (33). HHAT is an enzyme in the endoplasmic reticulum, which can palmitate hedgehog protein, and HHAT participates in the expression of the sonic hedgehog signaling pathway (34). Hedgehog is involved in almost all planned cell growth and survival mode of vertebrates, and plays an important role in the development of the forebrain and midbrain (35). Desmocollin 3 (DSC3) is a member of the desmocollin family (36), and DSC2 is a transmembrane protein that connects with the other three DSCs to ensure the integrity of myocardial structure and function and is also involved in cell signal transmission and apoptosis (37). The current research on DNA replication helicase/nuclease 2 (DNA2) is mainly focused on tumors, and there is currently no research on VaD. DNA2 is a mediator of genome stability and is required for proper Okazaki fragment maturation, for DSB repair, and for protection of replication forks. In addition, loss or mutation of DNA2 may inhibit DNA repair pathways in cells (38). Lipid metabolism plays an important role in vascular diseases such as atherosclerosis (39). A study has shown that high levels of receptor class A member 3 (SCARA3) methylation are associated with metabolic disorders, including type 2 diabetes and atherosclerosis, and SCARA3 may contribute to obesity and obesity related metabolic complications in an age-related manner (40). The potential role of these mRNAs in VaD needs to be further explored.

We also performed a GSEA analysis to determine the key regulatory pathways of lncRNA NEAT1. Natural killer (NK) cells are a type of cytotoxic lymphocyte that can infiltrate the ischemic brain and exacerbate neuronal cell death (41). Recent research indicates that NK cells play a key role in neurogenesis and cognitive decline during normal brain aging (42). It can be seen that lncRNAs may mediate the cytotoxicity of NK cells to aggravate the progression of VaD, and at the same time, the effect of NKT cells may be inhibited. Toll-like receptor (TLR) is an important pattern recognition receptor of the innate immune system, which triggers the inflammatory cascade by recognizing pathogens and related molecular patterns of damage (43). Toll-like receptor 2 (TLR2) and TLR4 are the two most important members of the TLR family. Previous studies have shown that TLR2 and TLR4 play a key role in inflammation after ischemic brain injury (44). Generally, TLR2 and TLR4 can be activated by a large number of exogenous and endogenous danger signal molecules, and ultimately promote the production of inflammatory cytokines, such as interleukin-6 (IL-6) and tumor necrosis factor-α (TNF-α) (45). Studies have shown that the TLR4 pathway induces the production of inflammatory cytokines, which leads to the development of cerebrovascular diseases during VaD (46). In addition, another study also showed that selective destruction of the TLR2-MyD88 interaction can inhibit inflammation in mice and improve their memory and learning ability (47). Therefore, reducing neuroinflammation by inhibiting the toll-like pathway may be an effective strategy for VaD.

In our research, for the first time, a ceRNA network related to a VaD immune response was constructed, and a sub-network related to NKT cells was constructed. However, our research still has certain limitations. We have not carried out any verification and mechanism experiments. In addition, there has been little research on the role of many of the mRNAs involved in VaD.

Conclusions

In conclusion, our research innovatively connects ceRNA to immunomodulation in the development of VaD. A decrease of NKT cells, and an increased expression of B4GALT1, PPP1R3B, MICB, HHAT, DSC2, DNA2, SCARA3, and lncRNA NEAT1 may have a significant correlation with the progression of VaD. These results will assist in a better understanding of the immunomodulatory mechanism of VaD from the perspective of ceRNA and may provide candidate therapeutic targets.

The article’s supplementary files as