Bioinformatics Analysis and Identification of Key Candidate Genes Associated with Nephrotoxicity Induced by Cobalt

Cobalt (Co) is a silver-gray, high-intensity, widely distributed metal element that exists in cobalt compounds, and its common valences are bivalence (Co²⁺) and trivalence (Co³⁺) ^[1]. The main routes of Co-exposure are occupational and environmental exposures. The human body can be exposed to high concentrations of Co²⁺ through inhalation of contaminated air, consumption of contaminated food and water, or ingestion of Co-containing supplements^[2], which has detrimental health effects when blood Co concentration reaches > 300 µg/L^[1]. The kidney is an important target organ, and animal experiments have shown^[3] that Co resulted in atrophy and degeneration of renal tubules and a significant upregulation of the kidney injury molecule 1 (Kim-1). In the presence of H₂O₂, Co²⁺ can participate in the Fenton reaction, generating highly reactive hydroxyl radicals that can bind to DNA, proteins, and lipids to promote tissue damage. However, the mechanism underlying Co nephrotoxicity is more complex and has not been clarified, creating the need for an in-depth and systematic study. With the advent of data sharing in recent years, data mining and the use of gene chips have become effective and necessary methods for toxicological studies. Previous studies on the mechanism of Co nephrotoxicity lacked systematic analysis involving the integration of bioinformatics and the identification of key pathways. A comprehensive and systematic study of gene expression profiles can help us better clarify the potential molecular mechanisms of Co nephrotoxicity.

The GSE57811 dataset used in this article is from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), with species of Rattus norvegicus. This dataset contains results on the renal effects of many toxicants. In this study, Co nephrotoxicity was investigated using a dataset consisting of three groups that were injected intraperitoneally with 10 mg/kg cobalt chloride (CoCl₂) every day continuously and collected on the fifth day. For the control group, to reduce the error, we selected all the data of the same batch, which were injected with an equal amount of saline intraperitoneally every day and collected on the fifth day for a total of 5 groups. The samples selected are shown in Supplementary Table S1 (available in www.besjournal.com).

The control samples in the GSE57811 dataset were used as the negative control group, and the exposed samples were used as the experimental group. Analysis of variance was performed using the limma package (version 3.44.3) in R language (version 4.0.2), the Bayesian method was used to build linear models, t-test and F-test were used to make comparisons between groups and obtain P-values, Benjamini-Hochberg correction was used for False Discovery Rate values and P-values, and multiples differences were obtained by log transformation. The differentially expressed genes (DEGs) of each group were filtered with a significance threshold of |logFC| > 1 and P < 0.05. After analyzing the GSE57811 dataset, 165 DEGs were identified, including 83 upregulated genes and 82 downregulated genes. The pheatmap (version 1.0.12) and ggplot2 (version 3.3.2) packages in R language were used to analyze and visualize the DEGs. Figures 1A and B show the volcano map and clustering heat map of DEGs, respectively.

Figure 1.  The results of bioinformatics analysis. (A) The volcano map of DEGs. The blue dots and red dots represent downregulated and upregulated genes, respectively; (B) The clustering heat map of DEGs;(C) GO enrichment analysis of DEGs; (D) KEGG pathway analysis of DEGs; (E) PPI network map of DEGs was constructed using the STRING online database and drawn using Cytoscape software; (F) PPI network map of top 5 hub genes identified by Cytohubba plug-in of Cytoscape, and their importance was represented by their color’s shade.

The obtained DEGs were uploaded to DAVID 2021 (https://david.ncifcrf.gov/home.jsp) for GO enrichment analysis and KEGG pathway analysis. Fisher’s exact test P-values were obtained, and a P-value < 0.05 indicated significance. GO enrichment analysis included three annotation contents: biological processes, molecular functions, and cellular components. The most significant GO-enriched biological processes of the DEGs included elastic fiber assembly, circadian rhythm, and response to drugs, as shown in Figure 1C. KEGG pathway analysis showed that DEGs were involved in circadian rhythm, steroid hormone biosynthesis, and FoxO signaling pathway, as shown in Figure 1D.

The PPI network graph was constructed using STRING 11.5 (https://cn.string-db.org/), setting the confidence level to 0.400 for screening. The results contained 143 nodes and 199 edges. The PPI network was visualized using Cytoscape (version 3.9.1) software, as shown in Figure 1E. The topological analysis was performed using the Cytohubba (version 0.1) plug-in of Cytoscape to screen the hub genes, among which the top 5 DEGs with the highest degree were selected as hub genes involved in Co-induced kidney injury. Based on the degree value, the top 5 hub genes were Basic Helix-Loop-Helix ARNT Like 1 [Arntl (Bmal1)], Neuronal PAS Domain Protein 2 (Npas2), Cyclin D1(Ccnd1), Nuclear Receptor Subfamily 1 Group D Member 1 (Nr1d1), and C-C Motif Chemokine Ligand 2 (Ccl2). The ranking by Degree method was shown in Supplementary Table S2 (available in www.besjournal.com), and their interaction relationships were shown in Figure 1F. Arntl, also known as Bmal1, is a basic helix-loop-helix protein, which is often considered to be part of the positive core clock. Npas2 is a transcriptional activator, and Npas2-Bmal1 can participate in the circadian feedback cycle by activating the transcription of biological genes. Ccnd1, an important member of the cell cycle protein family, is a positive cell cycle regulator, which can regulate cell cycle progression. Nr1d1, a member of nuclear receptor subfamily 1, negatively regulates the expression of core clock proteins, particularly for Bmal1^[4]. Ccl2 is a member of the CC subfamily of the chemokine family among inflammatory cytokines, which is secreted in proximal tubular epithelial cells. The most representative top 1 hub gene, Arntl, namely Bmal1, was selected to be validated via an animal experiment.

Though the species of the GSE57811 dataset was Rattus norvegicus, the toxicity mechanism is universal, and there will be no differences according to species. Therefore, ICR mice were selected for the verification experiments in this study. Thirty SPF-grade male healthy adult ICR mice weighing 20 ± 2 g were purchased from Changchun Yisi Experimental Animal Technology Co., Ltd. Mice were placed in the SPF-grade animal housing room with adequate water and food. Ethical approvals for all experimental protocols were obtained from the Ethics Committee of the School of Public Health at Jilin University.

Thirty mice, with adaptive feeding for 7 days, were randomly divided into 3 groups: saline control group, 16 mg/kg, and 20 mg/kg exposure groups (n = 10 each). Different concentrations of cobalt chloride (CoCl₂) solutions were intraperitoneally injected into the mice in the exposure groups, while the same volume of saline was intraperitoneally injected into the mice in the control group, and the mice in each group were injected daily for 14 consecutive days.

After 14 days of Co-exposure in each group of mice, the mice were humanely executed, and the kidneys were removed, fixed in 10% paraformaldehyde, paraffin-embedded, and then sectioned. The sections were dewaxed, rehydrated, and stained with hematoxylin and eosin (HE), and the pathological changes in the kidneys were observed under a microscope. The renal pathohistological results showed that the glomerular structure was normal in the control group, while the renal capsules were enlarged, the morphology of glomeruli was altered, and atrophy was obvious and accompanied by congestion in the 16 mg/kg and 20 mg/kg Co-exposure groups. In the 16 mg/kg and 20 mg/kg Co-exposure groups, edema was observed in the proximal convoluted tubules, and giant cells were found in the proximal convoluted tubule epithelium. The dilated tubules of the 20 mg/kg Co-exposure group displayed a hyaline cast, as shown in Figure 2A–C.

Figure 2.  Pathohistological observation of kidneys by HE staining and expression of Bmal1 protein in the kidneys by immunohistochemistry staining in mice exposed to different concentrations of CoCl₂. HE staining: (A) control group; (B) 16 mg/kg Co-exposure group; (C) 20 mg/kg Co-exposure group. The blue arrow points to the normal glomerular structure. The red arrows point to the enlarged renal capsules and the obvious shrunken glomeruli accompanied by congestion. The yellow arrows point to the proximal convoluted tubule edema and the proximal convoluted tubule epithelium giant cells. The black arrow points to the hyaline cast in the dilated tubules. Scale bar = 50 μm. Immunohistochemical detection of Bmal1 protein expression: (D) control group; (E) 16 mg/kg Co-exposure group; (F) 20 mg/kg Co-exposure group. The purple arrows point to immunohistochemical visualization of Bmal1 protein expression, and the color is brown. Scale bar = 50 μm. (G) Statistical graph of Bmal1 protein expression in the kidneys of mice exposed to different concentrations of CoCl₂ by immunohistochemistry. ^#P < 0.05, compared with the control group; ^*P < 0.05, compared with the 16 mg/kg Co-exposure group.

Kidney paraffin sections were dewaxed, rehydrated, and then subjected to antigen repair in an antigen repair solution at 98 °C for 5 min, followed by cooling to room temperature. The paraffin sections were incubated with 50 μL of endogenous peroxidase blocker (UltraSensitive^TM SP, KIT-9710, 2010279710G, MXB, Fuzhou, China) for 10 min and rinsed with PBS. Then, the sections were blocked with a non-specific staining blocker for 10 min, incubated with 50 μL of Bmal1 (14268-1-AP, ProteinTech, Wuhan, China, dilution concentration 1:300) per section overnight at 4 °C, and rinsed with PBS; subsequently incubated with the biotin-labeled secondary antibody dropwise for 10 min and rinsed with PBS; and then incubated again with streptavidin-peroxidase dropwise for 10 min, followed by rinsing with PBS. After DAB (DAB-0031, 2107120031E, MXB, Fuzhou, China) development and hematoxylin staining, photographs were taken. The expression of the Bmal1 protein was analyzed using Image J software (Version 1.53t), and the average optical density (AOD) indicated the protein expression. The data were analyzed using SPSS 24.0. The K-S test was used to test for normality; the data results following normal distribution were expressed as $\bar {\rm{x}}$ ± s. The means between multiple groups were compared using one-way ANOVA analysis, and statistical differences between each group were compared using Tukey’s method, both with a test level of α = 0.05. P < 0.05 was considered a statistically significant difference. The expression of Bmal1 protein in the kidneys of mice in the 16 mg/kg and 20 mg/kg Co-exposure groups was significantly decreased (P < 0.05) compared with the control group. The 20 mg/kg Co-exposure group had a lower protein expression of Bmal1 than the 16 mg/kg Co-exposure group (P < 0.05), as shown in Figure 2D–G. The immunohistochemistry results showed that the expression of Bmal1 protein in the kidneys of mice gradually decreased with the increase of Co dose, which was consistent with the bioinformatics analysis results.

Arntl, also known as Bmal1, is an important biological clock regulatory gene. Circadian disorders may be caused by the reduced expression of Bmal1^[5]. Studies by Nikolaeva et al.^[6] and Tokonami et al.^[7] have demonstrated that Bmal1 plays a role in maintaining normal renal functions. The loss of Bmal1 can result in impaired kidney and mitochondrial function^[6,8]. The renin-angiotensin-aldosterone system (RAAS) is influenced by circadian rhythms^[9], and a dysregulated RAAS increases the levels of angiotensin II (Ang-II), which activates the NF-kB pathway through the angiotensin-1 (AT₁) receptor and promotes the secretion of pro-inflammatory cytokines to induce renal inflammation^[5]. It was found^[7] that Bmal1^lox/lox/Ren1^dCre mice with inactivation of Bmal1 showed increasing urine output, a circadian rhythm disturbance of urinary Na⁺ excretion, and decreasing aldosterone levels in mice plasma. In addition, in the conditional knockout of Bmal1 (Bmal1^lox/lox/Pax8-rtTA/LC1) mice model, analysis of metabolomic and transcriptomic data and validation indicated that the biological clock in renal tubular cells controlled plasma urea nitrogen, arginine, and creatinine levels^[6]. The above research findings further support our speculation and verification in the animal experiment.

Based on the above results, we believe that one of the most important mechanisms of Co-induced renal injury is a change in circadian rhythm-related genes, especially the circadian rhythm gene Bmal1, which provides a new research target for Co-induced nephrotoxicity therapy. However, it still needs to be further explored by population and in vitro cellular experiments.