We used the keyword “diabetic heart” to search the microarray datasets from the GEO database of NCBI (https://www.ncbi.nlm.nih.gov/) (26). Ten pieces of datasets were displayed by restricting the entry type (datasets) and study type (expression profiling by array). We finally selected GSE4745 and GSE6880 as the microarray datasets by further filtration with the disease model (diabetic cardiomyopathy, DCM), sample (heart ventricles), title, summary and the purpose of this study. The GSE4745 dataset was based on the GPL85 platform, [RG_U34A] Affymetrix Rat Genome U34 Array done by Gerber et al. It contains 24 samples of the ventricular gene array data of rats (Rattus norvegicus), including 12 controls and 12 diabetics induced by streptozotocin (STZ) injection. In order to sufficiently observe the gene differences between control and diabetic group with typical DCM features, we selected eight samples including four controls and four diabetics on day 42 as the first set of microarray data for the analyses. The cardiac function in diabetic rats on that day showed significant ventricular diastolic and systolic dysfunction. The second dataset GSE6880 was based on the GPL341 platform, [RAE230A] Affymetrix Rat Expression 230A Array, including six samples of the ventricular microarray data of rats (Rattus norvegicus), of which, data from three normal hearts and three diabetic heart ventricles were selected as the second dataset for the analyses.

GSE4745 and GSE6880 microarray datasets were downloaded by using the GEOquery package (Version 2.54.1) of R software (https://www.rstudio.com/, Version 3.6.0); and were preprocessed before the screening of DEGs. Briefly, the data with null names or with multiple gene names were not used. As to that the same gene has multiple gene expression results, only the mean values were acquired. These collected data were normalized by utilizing “normalizeBetweenArrays” function, a quantile normalization in R, and the distributions for each set of microarray data were displayed by drawing a boxplot. Subsequently, the normalized microarray data were analyzed by using the limma package of R for DEGs identification. The limma package (Version 3.42.2) requires three matrices for differential analysis: expression matrix, design matrix, and contrast matrix. The design matrix was based on the annotation file of the microarray data and the experimental scheme. The contrast matrix was constructed according to the study groups. We set a threshold at logFC ≥1 (upregulated genes) or logFC ≤ −1 (downregulated genes) to screen out differentially expressed genes. The adjusted P < 0.05 was regarded statistically different. The results were saved in CSV file. The DEGs between normal and diabetic groups in datasets GSE4745 and GSE6880 were illustrated by Venn diagrams, from which the overlapped common up- and down-regulated DEGs can be obtained. The distribution and the cluster of the DEGs were visualized by volcano plots and heatmaps, respectively. The heatmap package (version = 1.0.12) was used to numerically normalize the expression levels of DEGs in GSE4745 and GSE6880 datasets in the direction of row normalization.

The GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses were performed by using the “clusterProfiler” package (Version 3.14.3) of the R language based on the common differentially up- and down-regulated genes of datasets GSE4745 and GSE6880, and an adjusted P < 0.05 was considered statistically significant. GO is an international standard classification system for gene functions, including biological process (BP), cellular component (CC) and molecular function (MF) (http://geneontology.org/) (27). KEGG is a database resource that integrates large scale of genomic, chemical and systemic functional information, including but not limited to molecular interaction and metabolic pathways, cellular processes, organismal systems related to human diseases and drug development (https://www.kegg.jp/).

Gene Set Enrichment Analysis (GSEA) was performed through the GSEA Pre-ranked commands in the GSEA software (http://www.broadinstitute.org/gsea/). The value of log2FC calculated by the limma package was used as the ranking metric. We used the C5 collection which contains gene sets annotated by GO terms, and the C2 pathway gene set collected from KEGG in the analysis. All genes from GSE4745 or GSE6880 were analyzed by GSEA separately, then we analyzed the overlapping pre-ranked enrichment terms which were both in top N (top50 for GO and top20 for KEGG pathway) terms of GSE4745 and GSE6880.

STRING database (https://string-db.org/, Version 11) is a worldwide online resource for the analysis of the interactions among known and predicted proteins. The PPI networks are viable tools to identify key gene modules and core hub genes between patients and healthy individuals. We constructed the PPI networks by importing the common DEGs into the STRING database. The interaction relationships among proteins encoded by these DEGs were then excavated between normal and DCM groups. Subsequently, we applied Cytoscape software to visualize the common DEGs on the basis of the PPI associations, and used the “cytoHubba” plug-in to obtain the hub genes. Top five hub genes were selected according to the ranking order of connectivity degree.

All studies were performed in accordance with the relevant guidelines and was approved by the Hebei University Animal Care and Use Committee (Approval No. IACUC-2019009XG). DCM model was established by using streptozotocin (STZ) treatment to the mice. Briefly, 9-week-old male C57BL/6N mice were subjected to intraperitoneal injections of STZ (Merck-Sigma-Aldrich®, USA, 100 mg/kg/day) dissolved in 0.1M sterile citrate buffer (adjusting the pH to 4.0 by using 1N NaOH) or vehicle for two consecutive days. Following 4-weeks, fasting blood glucose levels (overnight fasting for 12 hr starting around 9 pm to 9 am) were monitored using a glucometer (CONTOUR®PLUS Blood Glucose Monitoring System, Bayer, Germany). Mice with fasting blood glucose levels >11 mM were deemed diabetic followed by the echocardiographic analysis (28). The female mice were not used because estrogen may have an important role in glucose metabolism after STZ treatment and females are less sensitive to STZ than males (29–31).

Cardiac geometry and function were evaluated in anesthetized (isoflurane 1.2% administered with a calibrated vaporizer and an inhalant) mice using a two dimensional (2D) guided M-mode echocardiography (Vevo 2100) equipped with a 18–38 MHz linear transducer (Visualsonics, Canada). The heart was imaged in the 2D-mode in the parasternal long-axis view with a depth of 2 cm. The M-mode cursor was positioned perpendicular to interventricular septum and posterior wall of left-ventricle (LV) at the level of papillary muscles from the 2D-mode. Left ventricular diastolic and systolic anterior wall thickness (LVAWD, LVAWS), left ventricular diastolic and systolic posterior wall thickness (LVPWD, LVPWS), ejection fraction (EF), fractional shortening (FS), and LV mass were all measured. Normalized LV mass was calculated as LV mass (mg)/body weight (g).

At day 42 (week 6) following STZ injection, total RNA was isolated from ventricles using the TRNzol Universal Reagent (TIANGEN, China). RNAs were quantified using a NanoDrop^TM 2000c spectrophotometer (Thermo Fisher Scientific). Synthesis of cDNA and reverse transcription was performed using 1 μg total RNA in a 20-μl system following the instructions of FastKing RT Kit (With gDNase) (TIANGEN, China) and RT SuperMix for qPCR (APE × BIO, USA). Cyp1a1 and Rn18s (a housekeeping gene) primers for qPCR were designed by using the Primer5.0 software. Cyp2e1 and Col1a2 primers were referenced from the OriGene's website (https://www.origene.com.cn/). Col1a1, Col3a1 and Ctsk primers were referenced to the previous studies (32–34). The quantitative real-time PCR was performed by using a C1000 Touch Thermal Cycler CFX96^TM Real-Time System (BioRad) per the iQ^TM SYBR® Green Supermix (CLJC-Bio, China) instructions. Real-time PCR was triplicated for each cDNA sample. The primer sequences are shown in Table 1.

The procedures of mitochondrial protein extraction and Western blot was referenced to the previous published methods (35). Briefly, the heart ventricles were homogenized in an ice-cold MSE buffer containing 220 mM mannitol, 2 mM EGTA, 70 mM sucrose, 5 mM Mops (pH 7.4), 0.2% BSA and a protease inhibitor cocktail. The mitochondria pellet was obtained by a gradient centrifugation. The soluble mitochondrial protein fraction was produced by dissolving the mitochondria pellet in the lysis buffer [20 mM Tris/HCl (pH 7.4), 150 mM NaCl, 1 mM EDTA, 1 mM EGTA, 1% Triton, 0.1% SDS and 1% protease inhibitor cocktail] and centrifuged at 10,000 g for 30 min at 4°C (35). A portion (25 μg) of the mitochondrial protein was separated on 12% SDS-polyacrylamide gel, transferred electrophoretically to PVDF membranes and blotted against CYP4501A1 (1:1000, proteintech13241-1-AP), CYP4502E1 (1:4000, proteintech19937-1-AP) and GAPDH (loading control, 1:1000, Cell Signaling Technology) antibodies. After washing with TBST for three times, blots were incubated with horseradish peroxidase (HRP)-conjugated sheep anti-rabbit secondary antibody (1:5000, Cell Signaling Technology) at room temperature with a shaker for 1 h. Antigens were detected by the chemiluminescence method using Bio-Rad imager, and semi-quantitative band analysis was analyzed by Image Lab software (Bio-Rad, Version 5.1, ChemiDoc XRS) (4, 36).

The bioinformatic statistical analyses were performed in R software (Version 3.6.0). The codes have been uploaded to GitHub at https://github.com/ylchen0622/Paper-data-of-diabetic-cardiomyopathy-.git. The data based on animal experiments were presented as mean ± SEM. Statistical significance (p < 0.05) for each variable was estimated by an unpaired t-test (two-tailed). Softwares Origin was used for the correlation analysis between the log2-fold change of Cyp1a1 and Cyp2e1 mRNA levels and each of the cardiac function variables, as well as the correlation between log2-fold change of Cyp1a1 and fasting glucose level.

The distribution profiles of the microarray data in each dataset were matched via quantile normalization and boxplot analysis (Supplementary Figures 2A,B). A total of 212 DEGs (105 up-regulated and 107 down-regulated) in GSE4745, and 396 DEGs in GSE6880 (205 up-regulated and 191 down-regulated) were screened, among which, 37 common DEGs including 20 upregulated genes and 17 downregulated genes were identified in diabetic ventricles compared to normal controls (Figure 1). The names of these genes were listed in Table 2. The DEGs expression of each sample in the two datasets were visualized in the form of a heatmap in Figure 2. In addition, the results of microarray expression matrix in each dataset were displayed in the volcano plots (Figure 3). We screened three significantly down-regulated genes (LOC102553868, Hspa1a, Dbp) and eight significantly up-regulated genes (Acot1, Cyp26b1, Tgm1, Hmgcs2, Pdk4, Cyp2e1, S100a9, S100a8) in GSE4745 (Figure 3A); six significantly down-regulated genes (Card9, Adra1d, Crybb1, Sqle, Kazald1, LOC100909684) and three significantly up-regulated genes (Hmgcs2, Gal, Acot1) were screened in GSE6880 (Figure 3B). All comparison data between diabetic group and control (normal) group (diabetic vs control/normal) in datasets GSE4745 and GSE6880 were shown in the Supplementary Excel Files including upregulated and downregulated genes.

The GO functional enrichment analysis showed that the common up-regulated DEGs were mainly enriched in fatty acid metabolism during the biological process. At the cellular component level, they were mostly enriched on the inner mitochondrial membranes. In terms of the molecular function, they were chiefly enriched in the iron ions binding and heme binding functions (Figure 4A). Meanwhile, the analysis of KEGG metabolic pathway enrichment showed that these DEGs were mainly enriched in the chemical carcinogenic metabolic process and the metabolism of xenobiotics by cytochrome P450 (Figure 4B). On the contrary, the GO results revealed that common down-regulated DEGs were abundantly present in collagen fiber tissues, extracellular matrix tissues, and extracellular structure tissues. At the cellular component level, they were mostly amassed in collagen fiber trimers and extracellular matrix containing collagens. In respect to the molecular function, they were primarily enriched for the binding of platelet-derived growth factors and the binding of proteases and collagens (Figure 4C). At the KEGG metabolic pathway level, it is mainly enriched in the AGE-RAGE signaling pathway in protein digestion and absorption, and diabetes complications (Figure 4D). Along the same line, GSEA results showed that fatty acid metabolism, xenobiotic metabolism, and ATP synthesis through coupled electron transport chain in mitochondria are the main biological processes enriched in the overlapping terms of GSE4745 and GSE6880. The related pathways (p < 0.05) and involved genes, as well as the GSEA overlapping results including biological processes, molecular functions and cellular components were shown in the excel files as our Supplementary Materials.

We constructed a PPI network by introducing the screened common DEGs into the STRING online analysis software. By removing the unconnected nodes, there were a total of 37 nodes left in the network diagram (PPI enrichment p < 1.0^e−16), suggesting there are multiple biological interactions among these DEGs (Figure 5A). Subsequently, the results of the PPI analysis were imported into the Cytoscape software (Version 3.8.2) to screen out the hub genes. The first five hub genes were selected by using the cytoHubba plug-in function. The sequentially orders are as follows: Cyp1a1, Cyp2e1, Col1a1, Col3a1, Col1a2 (Figure 5B). These hub genes are closely connected and located at the hub of the PPI network, which are expected to become potential targeted gene candidates for the treatment of DCM.

Following 4 weeks of STZ treatment, the fasting blood glucose level reached to 275.4 mg/dL, which was nearly 2.3 folds higher than in control mice. The water and food consumption were comparably increased (Supplementary Table 1). STZ-induced diabetes caused significant reductions in the ejection fraction (EF), fractional shortening (FS), LV anterior wall thickness in systole (LVAWS) and LV posterior wall thickness in systole (LVPWS) (Figures 6A–D). Meanwhile, the LV end systolic diameter (LVESD) and left ventricular volumes in systole (LVVS) were significantly higher in mice with diabetes (Figures 6E,F). However, the stroke volume (SV), cardiac output (CO), the normalized LV mass, heart rate and related LV parameters in diastole including LVAWD, LVPWD, LVEDD and LVVD were not markedly altered following STZ injections (Figures 6G–N). The representative M mode and B mode images were shown in Supplementary Figure 3. Following 6 weeks of STZ treatment, the liver mass and size were significantly increased, while the spleen mass and size were dramatically reduced compared to the control. STZ failed to affect the body weight, heart weight and kidney weight themselves. However, the kidney to body weight was shown to be increased (Supplementary Table 2).

Levels of hub genes including Cyp1a1, Cyp2e1, Col1a1, Col3a1 and Col1a2, as well as Ctsk were measured in heart ventricles of control and STZ induced diabetic mice using quantitative real-time PCR. Cyp1a1 and Ctsk mRNA expressions were elevated by about 12-fold change and 2-fold change, respectively. However, Cyp2e1, Col1a1, Col3a1 and Col1a2 had no obvious difference in diabetic mice than non-diabetic subjects (Figures 7A–F).

Our analyses showed that Cyp1a1 had a positive correlation with ejection fraction, fractional shortening, and LVPWS in non-diabetic mice. On the contrary, a negative correlation was illustrated between them in diabetic mice (Figures 8A–C). At the same time, it revealed a negative correlation between Cyp1a1 and LVESD or LVVS in non-diabetic mice but a positive correlation in diabetic mice (Figures 8D,E). Additionally, we observed negative correlations between Cyp1a1 and LVAWs in both diabetic and non-diabetic groups (Figure 8F). Similarly, there was a positive correlation between Cyp2e1 and EF, FS or LVPWS in non-diabetic mice, while a negative correlation between them in diabetic mice (Figures 8G–I). Cyp2e1 was negatively correlated with LVESD and LVVS in non-diabetic mice but positively correlated with them in diabetic mice (Figures 8J,K). However, we observed positive correlations in both groups between Cyp2e1 and LVAWs (Figure 8L). In addition, correlations between Cyp1a1 mRNA levels and fasting glucose in non-diabetic and diabetic mice were also shown in Supplementary Figure 4.

Protein levels of mitochondrial CYP4501A1 and CYP4502E1 were measured in heart ventricles of non-diabetic and STZ induced diabetic mice by using immunoblot analysis. CYP4501A1 expression in mitochondria were significantly increased in diabetic heart compared with control. While there was no significant change of mitochondrial CYP4502E1 expression between two groups although there was an increased trend induced by STZ (Supplementary Figure 5).