Identification of key genes and therapeutic drugs for cocaine addiction using integrated bioinformatics analysis

Data

Cocaine addiction-related data GSE54839 (Homo sapiens), GSE67281 (Homo sapiens), GSE186981 (Mus musculus), and GSE155313 (Mus musculus) were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). The dataset GSE54839 was chosen as our experimental set, which is based on the GPL6947 platform (Illumina HumanHT-12 V3.0 expression beadchip). This microarray-based study determined the profiles of midbrain gene expression in chronic cocaine abusers (n = 10) and well-matched drug-free control subjects (n = 10). Array-related procedures were performed in triplicate for each subject.

GSE67281 (Homo sapiens), GSE186981 (Mus musculus), and GSE155313 (Mus musculus) were chosen as our validation sets. GSE67281 is an expression profile in postmortem human midbrain specimens from chronic cocaine abusers (n = 11) and well-matched control subjects (n = 11). The GSE186981 is RNA-Seq data in hybrid mouse diversity panel (HMDP) mouse strains of nucleus accumbens (NAc) and prefrontal cortex brain regions (PFC). GSE155313 is the RNA-Seq data from the VTA region of mouse that underwent one of four commonly used paradigms: acute home cage injections of cocaine, chronic home cage injections of cocaine, cocaine-conditioning, or intravenous-self administration of cocaine.

Human PPI data were downloaded from the STRING database (https://string-db.org/) (Szklarczyk et al., 2021). With a combined score of >900 as the threshold, a total of 230,524 interactions between 11,763 genes were obtained.

Cocaine addiction-related PPI network

The dataset GSE54839 was differentially analyzed using the R package “limma” to obtain their differential genes. The p-value of <0.05 and |log₂F>p20mm| > 0.2630344 (i.e., fold change ≥ 1.2 or fold change ≤ 0.8) were considered statistically significant.

To obtain the interaction relationships between DEGs, the downloaded PPI data were filtered using the DEGs obtained from the microarray data GSE54839. By using the gene interactions as edges and DEGs as nodes, a differential gene network was constructed. After removing scatters from the network, the core network was defined as a PPI network related to cocaine addiction.

Key gene identification

We proposed a centrality algorithm integration strategy to analyze genes in cocaine addiction-related PPI networks. The scores of the node in network under each centrality algorithm were calculated separately by applying a series of centrality measures, including degree, edge-percolated component (EPC), Laplacian centrality, maximum neighborhood component (MNC), Katz radiality, and semi-local centrality (SLC). The intersection of the top 10 genes of each centrality algorithm was considered the key gene.

In this study, G is the cocaine addiction-related PPI network we built, and V(G) is the collection of nodes in the network. For node x in G, N(x) is the set of direct neighbors of x in G. For collection A, |A| is used to represent the number of elements in the collection. The specific algorithms are as follows:

1. Degree (Deg)

   ()=|()|.
2. Edge percolated component (EPC)

   ()=1|()|∑=11000∑∈.

   Given a threshold of 0.5, 1,000 reduced networks were created by assigning each edge a random number between 0 and 1 and removing edges with associated random numbers less than the threshold. Let the G_k be the reduced network generated at the kth reduced process. If nodes x and y are connected in G_k, set  to 1; otherwise, =0 (Chin et al., 2014).

3. Laplacian centrality (Qi et al., 2012):

   ()=()2+()+2∑∈()().

   (1)
4. Maximum neighborhood component (MNC):

   ()=|((())|,

   (2)

   where m(N(x)) is a maximum connected component of the induced subgraph of G by N(x) (Lin et al., 2008).

5. Katz centrality:

   ()=∑=0∞∑=1|()|(),

   (3)

   where A is the adjacency matrix of the network G with eigenvalues λ, () is the number of paths from x to y with length k, α is a damping factor and 0<<1max. In all our experiments, we chose α = 0.1 (Wei et al., 2020).

6. Radiality (rad):

   ()=∑∈()+1−(,)|()|−1,

   (4)

   where d is the diameter of the network G (Valente and Foreman, 1998).

7. Semi-local centrality (SLC):

   ()=∑∈()∑∈()(),

   (5)

   where B(z) is the number of direct connections and two-step neighbors for node z (Chen et al., 2012).

Correlation analysis of key genes with cocaine addiction

The Human Protein Atlas (HPA; https://www.proteinatlas.org/) database creates a brain-centric knowledge resource on RNA and protein expression in three mammalian brains: human, pig, and mouse (Sjöstedt et al., 2020). The RNA expression of key genes in different brain regions in humans and mice was searched in the brain section of the HPA database.

The Comparative Toxicogenomics Database (CTD, https://ctdbase.org/) was used to obtain associations between key genes and cocaine addiction. In the CTD database, the inference score reflects the degree of similarity between the CTD chemical–gene–disease network and a similar scale-free random network (Davis et al., 2023). The higher the score, the higher the degree of association between the disease and the gene.

The role of key genes in the mechanism of cocaine addiction was identified by text mining in the PubMed database. The search keywords were “cocaine addiction” and the four key genes.

Three independent sets GSE67281 (human), GSE186981 (mouse), and GSE155313 (mouse) were used for pre-addiction and post-addiction differential expression analyses to verify key genes, and the threshold and differential analysis methods were consistent with those of the experimental set GSE54839.

The R package “homologene” was used to search for homologous genes between the human and the mouse. The “homologene” package is a package based on the NCBI HomoloGene (https://www.ncbi.nlm.nih.gov/homologene/) database. The HomoloGene database is a system that can automatically detect congeners in human and mouse genes (NCBI Resource Coordinators, 2014).

To investigate the possible molecular mechanisms of key genes for cocaine addiction, we used the Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.kegg.jp/ or http://www.genome.jp/kegg/) database for enrichment analysis. An adjusted p-value of <0.05 was considered to be statistically significant.

Key gene

The scores of each node in the cocaine addiction related PPI network were calculated separately using seven different centrality algorithms (Figure 3A, Supplementary Table 2). Based on our proposed centrality algorithm integration strategy, the top ten genes scored by each algorithm were selected, and their intersections (FOS, IL6, EGR1, and JUN) were regarded as key genes (Figure 3B). The scores of the seven centrality algorithms for the four key genes are shown in Figure 3C.

We performed correlation analyses for four key genes in non-drug control participants (Figure 3D) and chronic cocaine abusers (Figure 3E), respectively. The correlation heat map showed that all four key genes were positively correlated, and the correlation showed a significant increase in the cocaine group. It suggested that these four key genes might be more closely related to each other and had synergistic effects during addiction. In addition, the expression of the IL6 gene was the lowest of the four key genes.

Expression of key genes in brain regions

To investigate the expression of four key genes in the brain, we searched the HPA database for the expression of four key genes in different brain regions in the human and mouse (Figures 4A, ​,B).B). The results showed that four key genes were expressed in all regions of the human brain, and IL6 was expressed in the human midbrain region lower than the other three genes, which is consistent with our findings. IL6 was not detected in the mouse midbrain, hypothalamus, pituitary gland, retina, pons, and medulla.

Correlation analysis of CTD databases

The correlation scores between each gene in the cocaine addiction-related PPI network and cocaine addiction were searched in CTD. The scores for the top 100 genes are shown in Figure 5. It showed a higher degree of association between the four key genes and cocaine addiction. CTD showed that FOS and EGR1 could be biomarkers of cocaine addiction or play a role in addiction, and EGR1 could be a gene for a therapeutic target in the treatment of cocaine addiction.

Literature validation

The PubMed database showed that four key genes were all associated with cocaine addiction. Both FOS (Fos proto-oncogene) and JUN (Jun proto-oncogene) are members of the AP-1 transcription factor complex. Multiple studies have shown that cocaine affected the expression of FOS proteins (Todtenkopf et al., 2002; Imam et al., 2005; Larson et al., 2010; Lobo et al., 2010). Zhang et al. (2004) prepared CPu extracts from D1 and D3 receptor mutant mice and wild-type control littermates at different time points after cocaine injection and found that ERK activation mediates acute cocaine-induced expression of c-fos (Fos). The study by Xu (2008) found that c-fos might mediate cocaine-induced persistent changes by regulating the formation of AP-1 transcriptional complexes and gene expression. Previous studies have demonstrated that cocaine causes increased expression of the JUN protein (Malaplate-Armand et al., 2005; Paletzki et al., 2008). Cocaine affects the expression of the JUN protein (Imam et al., 2005).

There are some studies proving that EGR1 (early growth active protein 1) and c-fos expressions are reduced after cocaine induction (Helton et al., 1993; Ennulat et al., 1994). In experiments on mutant mice by Valjent et al. (2006), EGR1 was found to play a vital role in cocaine-related behavior. Humblot et al. (1998) found that acute cocaine administration was effective in inducing c-FOS and EGR-1 direct early genes, and cocaine-induced EGR-1 and c-FOS expression was significantly reduced in brain regions of rats.

IL6 (interleukin-6) is a pro-inflammatory cytokine. The study by Halpern et al. (2003) showed that men and women respond weakly to pro-inflammatory challenges to IL6 after intravenous cocaine. In experiments measuring changes in IL6 levels in crack cocaine-dependent adolescents after 21 days of withdrawal, it was found that IL6 was elevated in patients on admission compared to the control group (Pianca et al., 2017).

Independent set analysis

A differential expression analysis was performed on the human dataset GSE67281, which is the expression profile of human cocaine abusers in the midbrain region. A total of 200 DEGs were identified after annotation, including 110 upregulated genes and 90 downregulated genes. There were 20 intersecting genes in the datasets GSE54839 and GSE67281 (Figure 6), including the key genes JUN, FOS, and EGR1, all of which were downregulated in the addictive state (Table 1), while the IL6 gene was not annotated in GSE67281.

A differential expression analysis and a homology analysis were performed on the mouse dataset GSE155313 from the VTA brain region, and the Fos gene was identified as a downregulated differential gene under four different conditions. The degree of difference in the Fos gene was not the same between chronic and acute home cage injections of cocaine, and even greater in chronic home cage injections of cocaine condition (Table 2). It suggested that the Fos gene plays a crucial role in long-term addiction.

The JUN, FOS, and EGR1 genes were shown to be downregulated differential genes in the human validation set, and the FOS gene was also downregulated in the mouse validation set, which is the same as the experimental set. Analysis of the independent sets showed that the key genes we identified were well-confirmed.

Pathway verification

To reveal the roles of the key genes, we performed a KEGG enrichment analysis of all the genes in the cocaine addiction-related PPI network. A total of 75 KEGG pathways were enriched (Figure 7A). The four key genes were mainly enriched in the TNF signaling pathway, cocaine addiction, amphetamine addiction, IL-17 signaling pathway, MAPK signaling pathway, and Toll-like receptor signaling pathway. Previous studies have shown that these pathways were all linked to cocaine addiction (Northcutt et al., 2015; Lewitus et al., 2016; Brown et al., 2018; Ganguly et al., 2019; Montesinos et al., 2020; Bingor et al., 2021). The subnetwork associated with the key genes and their enrichment pathways (Figure 7B) showed that the four key genes were closely connected in the network, among which JUN, FOS, and IL6 were enriched into multiple pathways, and EGR1 was closely related to these pathways and played a very important synergy.

The cocaine addiction pathway and the amphetamine addiction pathway are two enriched addiction-related pathways (Figures 8A, ​,B).B). They have similar addiction mechanisms, and both have enhanced firing activity of dopamine neurons in the VTA of the midbrain, resulting in enhanced dopamine release from the NAc.

With the stimulation of addictive drugs, the FOS gene induces and maintains an addictive state in the short term of addiction. In human who achieve long-term addiction after further drug use, the FOS gene causes long-term adaptive changes in the brain, and the JUN gene dimerizing with ΔFosB leads to an increased cocaine response. ΔFosB desensitizes c-fos mRNA induction after chronic amphetamine exposure (Renthal et al., 2008). Zhang et al. (2006) have shown that FOS might mediate cocaine-induced persistent changes by regulating AP-1 transcriptional complexes and target gene expression. To sum up, both our key genes JUN and FOS played important roles in the addiction pathways.