We obtained the GSE76894 and GSE76895 through the NCBI GEO public database (http://www.ncbi.nlm.nih.gov/geo/). Data for 19 T2DM patients and 84 control samples from GSE76894, 36 T2DM patients, and 32 control samples from GSE76895. Two microarray data were based on the GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix, Santa Clara, CA, United States).

Two OP-related datasets, GSE56815 and GSE7429, were also selected from the GEO website, which was done in the same way. The samples in the GSE56815 and GSE7429 datasets are osteoporosis patients characterized by low BMD. Dataset GSE56815 includes 40 female controls (20 premenopausal, 20 postmenopausal) and 40 female patients with low BMD (20 premenopausal, 20 postmenopausal). Dataset GSE7429 includes 10 postmenopausal controls (age 56.8 ± 1.93 years) and 10 postmenopausal patients (age 57.9 ± 1.45 years) with low BMD. Microarray data were based on the GPL96 [HG-U133A] Affymetrix Human Genome U133A Array (Affymetrix, Santa Clara, CA, United States).

We combined two datasets for each disease to increase the sample size. R software (version 4.2.3; https://www.r-project.org/) and used the “limma” package in the R software to screen for DEGs between the control and T2DM groups and between the control and OP groups, respectively, with parameters set to adj. p-value <0.05. We then used the “ggplot2” and “heat map” packages to draw volcano maps and heat maps. DEGs with log2FC > 0 were considered to be upregulated genes, while the downregulated genes were screened according to log2FC < 0.

WGCNA, a method for identifying highly co-expressed gene modules and their relationship with disease (Rasmussen et al., 2020), utilizes WGCNA software package R.4.2.3 for analysis. In this study, the WGCNA R package was used to construct coexpression networks integrating clinical characteristics for the DEGs of T2DM and OP. Initially, hierarchical cluster analysis using the Hclust function in R was conducted to eliminate outlier samples. Subsequently, an adjacency matrix was established using a soft-threshold β and a gene–gene correlation matrix to describe the degree of association between the nodes. The resultant clusters were analyzed via a topological overlap matrix, assigning modules by color and module feature (ME). Furthermore, Pearson’s correlation test was applied to calculate the correlation between ME and clinical features. The modules with a |ME|> 0.3 and a p-value <0.05 were deemed to have significant interactions with clinical features (Zhang et al., 2023).

The potential shared genes were identified as T2DM-related DEGs overlapped with the OP-related ones. These shared genes could have the potential ability to link the pathogeneses of T2DM and OP.

To further determine the biological features of potentially shared genes, Gene Ontology (GO) analysis was accomplished by the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool. The mission of the DAVID knowledgebase is to support biomedical and other research by providing comprehensive information about the evolution of protein-coding gene families, particularly protein phylogeny, function, and genetic variations impacting that function. GO enrichment analysis was performed using the DAVID online tool, including molecular functions (MFs), cellular components (CCs), and biological processes (BPs). The Encyclopedia of Genes and Genomes (KEGG) was used for pathway enrichment analyses (https://bioinfo.org/kobas/). The enriched function with p < 0.05 was considered a significant pathway. Based on this analysis method, we selected the top 15 GO biological processes and 10 KEGG pathways.

To better screen the risk of crosstalk genes between T2DM and OP, LASSO regression was performed in R project. The Co-DEGs of T2DM and OP were retained for feature selection, and the best tandem genes were initially identified using LASSO regression. To further narrow the scope, we combined the two results of this LASSO algorithm by taking the intersection set.

Six-week-old male C57BL/6J WT mice were randomly divided into two groups (n = 5 or 6 mice/group): (1) standard diet control group (SD) and (2) HFD + streptozotocin-treated group (HFD + STZ). Mice in the SD group were fed a standard chow diet. In the HFD + STZ group, mice were first fed an HFD (#D12492) for 4 weeks and then injected with STZ (30 mg/kg, Sigma-Aldrich) for three consecutive days to induce T2DM, as previously described (Ma et al., 2020). The mice with blood glucose levels over 16.7 mmol/L measured by a blood glucose meter (LifeScan) were considered to have T2DM and were continually an HFD until 16 weeks.

Eight-week-old female C57BL/6J WT mice were randomly divided into two groups (n = 5 or 6 mice/group): 1) sham surgery and 2) ovariectomy (OVX). The mice were anesthetized by intraperitoneal injection with 1% sodium pentobarbital solution (50 mg/kg), and the bilateral ovaries were removed through a dorsal approach (Shen et al., 2020). After finding the ovary, the ovaries of the mice in the sham group were returned to the abdominal cavity without resection.

All mice were housed in a specific pathogen-free (SPF) animal laboratory in a 12-h light:12-h dark cycle at room temperature (20°C–22°C). All experiments were performed according to procedures approved by the Animal Research Committee of Tongji Medical College, Huazhong University of Science and Technology, Hubei Province, China (Reference Number: 3,053, Date Approved: 10 October 2020).

Tail fasting blood glucose was assessed in mice using a blood glucose meter (LifeScan) following a 12-hour fast. For the intraperitoneal glucose tolerance test (IPGTT), mice underwent injection with glucose (2 g/kg body weight, i.p.) after a 12-hour fast, with blood samples collected from the tail vein at 0 min, 30 min, 60 min, 90 min, and 120 min post-injection. Similarly, the intraperitoneal insulin tolerance test (IPITT) involved injecting insulin (0.75 U/kg body weight, i.p.) into fasted mice, with blood samples collected from the tail vein at 0 min, 15 min, 30 min, 60 min, and 90 min thereafter.

Paraffin tissue paraffin sections were subjected to deparaffinization with xylene, followed by dehydration in a series of 100%, 95%, 85%, and 70% ethanol. The sections were incubated overnight at 4°C with antibodies specific for anti-insulin (1:100, Abcam, U.K.), glucagon (1:100, Abcam, U.K.), vanin1 (1:100, Proteintech, China), and osteopontin (1:100, R&D Systems, United States). After three washes with PBS for 10 min, the sections were incubated with fluorescent secondary antibody (1:100 dilution) for 1 h, and nuclei were counterstained with DAPI (Servicebio, Wuhan, China) at room temperature for 10 min in the dark (Chen et al., 2023a). The sections were observed by a fully automated slice scanning system (VS120, Olympus, Japan).

Total protein in cell lysates was extracted using a radio-immunoprecipitation assay (RIPA) buffer containing protease and phosphatase inhibitors. Tissue lysates were boiled in SDS loading buffer for 10 min; proteins were then resolved by SDS-PAGE and transferred to polyvinylidene difluoride membranes (PVDF, Millipore Corp, MA, United States). After incubation with 5% BSA in TBST for 1 h, the membranes were incubated with anti-VNN1 antibody (1:1000, Proteintech, China) at 4 °C overnight. The membranes were washed three times with TBST and then incubated with horseradish peroxidase-conjugated anti-rabbit antibody (1:2000 dilution) for 1 h.

The three-dimensional structure of trabecular bone was analyzed in the femur on a micro-CT scanner (Skyscan-1176, Bruker, Karlsruhe, Germany). The system was configured with a 1-mm AL filter, 9 μm resolution, 50 kV voltage, and 200 μA current. For the femur, scanning started at the lower growth plate and extended proximally for 600 slices, which included trabecular bone (slices 1–400) and cortical bone (slices 500–550). The following parameters were measured by CT analyzer (version 1.15.4.0): bone volume (BV), total volume (TV), trabecular bone volume fraction (BV/TV), trabecular bone thickness (Tb.Th), trabecular bone number (Tb.N), trabecular bone separation (Tb.Sp), structural model index (SMI), and cortical bone thickness (Cb.Th).

Total RNA was extracted from mouse pancreas and bone tissues using TRIzol reagent, and cDNA was reverse-transcribed using the HiScript III RT SuperMix (Vazyme, China) under the following conditions: 15 min at 37 °C and 5 s at 85 °C. Next, the cDNA was amplified using ChamQ SYBR qPCR Master Mix (Vazyme, China). β-Actin was used as an internal control in 10 μL volume. The 2^−ΔΔCq method was used to calculate the relative mRNA levels. The relative mRNA levels were calculated by the 2^−ΔΔCq method, and the primer sequences are shown in Table 1.

The VNN1 levels were determined using the Human Mouse Vanin 1 (VNN1) ELISA Kit (Jianglai, China) according to the manufacturer’s instructions. Absorbance at 450 nm was detected using an enzyme labeler (PerkinElmer, Waltham, MA, United States).

The clinical samples involved in the experiment were all from our previous studies (Chen et al., 2023b). The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Ethics Committee of Tongji Medical College, Huazhong University of Science and Technology (ChiCTR2000034751). All subjects gave their informed consent for inclusion before they participated in the study.

Type 2 diabetes mellitus group: 1) typical diabetes mellitus symptoms (irritable thirst and excessive drinking, polyuria, polyphagia, and unexplained weight loss); 2) random blood glucose ≥11.1 mmol/L; 3) fasting blood glucose ≥7.0 mmol/L; 4) blood glucose after glucose loading for 2 h ≥ 11.1 mmol/L; and 5) glycated hemoglobin ≥6.5%, with two of the above conditions fulfilled, and exclusion of other types of diabetes.

Based on the diagnostic criteria of the International Society for Clinical Densitometry (Watts et al., 2013), we measured BMD in participants by dual-energy X-ray absorptiometry (DXA) and calculated T-scores. We then categorized participants into normal (T-score ≥ −1), osteoporotic (−2.5 < T-score < −1), and osteopenia (T-score ≤ −2.5) according to T-scores. In this study, osteoporosis was diagnosed as femoral neck BMD ≤ −2.5 SD measured by DXA.

Type 2 diabetes mellitus combined with osteoporosis (T2DM + OP) group: 1) meet the diagnosis of type 2 diabetes mellitus; 2) meet the diagnosis of osteoporosis.

The data were analyzed using the GraphPad Prism software (8.0.1). Differences in numerical parameters between the two groups were assessed by an unpaired two-tailed t-test. All data are expressed as the mean ± SD, and p < 0.05 was considered statistically significant.

Two microarray datasets (GSE76894 and GSE76895) of T2DM were collected from the GEO database. The datasets comprised data from 19 T2DM patients and 84 control samples from GSE76894 and 36 T2DM patients and 32 control samples from GSE76895. To identify DEGs between T2DM and healthy controls, we gathered microarray expression profiles of GSE76894 and GSE76895 from the GEO database. Following the consolidation and normalization of the microarray data, we identified 8,506 DEGs between T2DM and healthy controls utilizing the “limma” package (p < 0.05), comprising 4,121 upregulated genes and 4,385 downregulated genes. Subsequently, we generated a volcano map (Figure 1A) and a hierarchically clustered heat map of the DEGs (Figure 1B).

We also selected two OP datasets (GSE56815 and GSE7429) from the GEO database and processed them with the same method, which resulted in 2,030 DEGs, including 1,091 upregulated genes and 939 downregulated genes. The results were plotted as a volcano map (Figure 1C) and a heat map (Figure 1D).

A WGCNA was used to detect clusters of co-expressed genes exhibiting expression differences between T2DM and OP and to assess the correlation of combined modules with disease phenotypes. The T2DM model had a soft-threshold β = 16, whereas the OP model had a soft-threshold β = 8 (Figures 2A, B). After combining similar gene modules, nine modules were recognized in the T2DM model set and 11 modules in the OP model set. As shown in Figure 2C, the brown module had the highest positive correlation with the occurrence of T2DM (r = 0.31), and the turquoise module had a substantial negative correlation with the occurrence of T2DM (r = −0.35). Furthermore, in the OP modeling set, the gray module demonstrated the strongest positive correlation with OP occurrence (r = 0.35), and the red module showcased the most prominent negative correlation with OP occurrence (r = −0.3) (Figure 2D).

A total of 19 genes overlapped in the strongest positive and negative modules of T2DM and OP. These genes may be related to the pathogenesis of T2DM and OP. We used the DAVID online tool to perform GO analysis of these 19 co-DEGs, including biological processes (BPs), cellular components (CCs), and molecular functions (MFs) (Figures 3A–C). The BPs were mainly focused on the cell–cell adhesion, positive regulation of I-κB kinase/NF-κB signaling, signal transduction, negative regulation of cell proliferation, and immune response; CCs were mainly focused on the integral component of the membrane and the Golgi, plasma, and lysosomal membranes; and MFs were mainly focused on zinc ion binding, RNA polymerase II core promoter proximal region sequence-specific DNA binding, metal ion binding, protein binding, and tumor necrosis factor receptor binding. We then performed KEGG enrichment analysis using the KOBAS website and identified 10 KEGG pathways (Figure 3D). The KEGG pathway analysis showed that the DEGs were closely related to pantothenate and CoA biosynthesis and the glycosaminoglycan biosynthesis–chondroitin sulfate/dermatan sulfate, renin–angiotensin system, and other pathways.

We then extracted the expression data of the 19 genes from the T2DM gene expression profile using the LASSO algorithms with an optimal lambda value of lambda.min = 0.0055 (Figures 4A, B). The results from the LASSO regression identified 12 genes as signature genes in T2DM. We used the same method to obtain eight genes for the diagnosis of OP (lambda.min = 0.0385; Figures 4C, D). The diagnostic genes of the two diseases were intersected to obtain six common diagnostic genes, namely, LGALS2, LTB, STEAP4, TPBG, VGF, and VNN1.

We further validated the changes of hub genes in islets of type 2 diabetic mice. To induce a typical T2DM model, a low dose of STZ was consecutively injected intraperitoneally for HFD-fed mice. The T2DM group had significantly increased body weight, and fasting plasma glucose was significantly higher than that of the SD group (Figures 5A, B), which showed obvious impaired glucose tolerance and insulin resistance (Figures 5C, D). In addition, mice in the T2DM group had an increased number of pancreatic islet α-cells and a disproportionate α/β-cell ratio (Figures 5E, F). Interestingly, PCR results from the mouse pancreas showed that LGALS2, LTB, TPBG, and VNN1 expression was upregulated, and VGF expression was downregulated in the T2DM group of mice compared with the SD group (Figure 5G).

We used the OVX mice model to simulate osteoporosis in postmenopausal women. Micro-CT was used to analyze the trabecular bone changes in the distal femur of different model groups 10 weeks after OVX. Compared with the sham group, the OVX  group demonstrated a significant decrease in BV/TV, Tb.Th, Tb.N, and Cb.Th and an increase in Tb. Sp and SMI (Figures 6B–G). PCR results derived from mouse femurs showed that LTB, STEAP4, TPBG, and VNN1 expression was upregulated in the OVX group of mice compared with the sham group (Figure 6H).

Through bioinformatics methods, we found that VNN1 was an upregulated DEG in islet samples of T2DM and OP. VNN1 is an emerging biomarker that has not yet been detected in diabetic osteoporosis patients, so we chose VNN1 for further study. Immunofluorescence results indicated that VNN1 expression was upregulated in the T2DM group (Figures 7A–C). Western blot showed the same results, revealing a substantial increase in VNN1 protein expression in T2DM groups compared to the SD group (Figures 7E, F). To clarify the distribution location of VNN1 in the mouse pancreas, we co-stained VNN1 with insulin (INS) and glucagon (GCG). We found that VNN1 was mainly distributed in pancreatic islets and partially merged with both INS and GCG. We also examined the expression of VNN1 in the bones of mice and found that the expression of VNN1 was significantly increased in the bones of T2DM model mice (Figures 7G–K). We used the osteogenic marker osteopontin (OPN) for immunofluorescence double staining with VNN1. Zooming in on the images, we found that VNN1 was localized around osteoblasts (red rectangles), bone marrow adipocytes (yellow circles), and blood vessels (white rectangles) in the bones of T2DM model mice (Figure 7H). Similarly, serum VNN1 expression was increased in T2DM mice compared to controls (Figure 7D).

Next, we investigated the expression of VNN1 in the pancreas and bone of OVX mice. Similarly, pancreatic VNN1 expression was increased in OVX mice (Figures 7E, F), localized predominantly in the pancreatic islets, and partially merged with INS and GCG (Figures 7A–C). VNN1 expression was increased in the serum of OVX mice compared with the sham group (Figure 8D). Immunofluorescence showed that VNN1 expression was increased in the bones of OVX mice, distributed around osteoblasts (red rectangles), bone marrow adipocytes (yellow circles), and blood vessels (white rectangles) (Figure 8H). Western blotting showed the same results, with a significant increase in VNN1 protein expression in the OVX group compared with the sham group (Figures 8J, K).

In order to verify the clinical application potential of the VNN1 gene, the ELISA method was used to detect the protein level encoded by the VNN1 gene in clinical samples. The T2DM group and the T2DM combined OP (T2DM + OP) group had higher blood glucose levels than the control group; the DXA femoral T-score was lower in the OP group and the T2DM + OP group than the control group, and the difference was statistically significant (Table 2). The expression of serum VNN1 in the T2DM, OP, and T2DM + OP groups was higher than that in the control group, and the difference was also statistically significant (Figure 9). These findings indicate that VNN1 is a potential biomarker for the early diagnosis and prognostic prediction of T2DM and OP.