Genome-wide differential expression profiling of lncRNAs and mRNAs associated with early diabetic cardiomyopathy

Diabetic cardiomyopathy is one of the main causes of heart failure and death in patients with diabetes. There are no effective approaches to preventing its development in the clinic. Long noncoding RNAs (lncRNA) are increasingly recognized as important molecular players in cardiovascular disease. Herein we investigated the profiling of cardiac lncRNA and mRNA expression in type 2 diabetic db/db mice with and without early diabetic cardiomyopathy. We found that db/db mice developed cardiac hypertrophy with normal cardiac function at 6 weeks of age but with a decreased diastolic function at 20 weeks of age. LncRNA and mRNA transcripts were remarkably different in 20-week-old db/db mouse hearts compared with both nondiabetic and diabetic controls. Overall 1479 lncRNA transcripts and 1109 mRNA transcripts were aberrantly expressed in 6- and 20-week-old db/db hearts compared with nondiabetic controls. The lncRNA-mRNA co-expression network analysis revealed that 5 deregulated lncRNAs having maximum connections with differentially expressed mRNAs were BC038927, G730013B05Rik, 2700054A10Rik, AK089884, and Daw1. Bioinformatics analysis revealed that these 5 lncRNAs are closely associated with membrane depolarization, action potential conduction, contraction of cardiac myocytes, and actin filament-based movement of cardiac cells. This study profiles differently expressed lncRNAs in type 2 mice with and without early diabetic cardiomyopathy and identifies BC038927, G730013B05Rik, 2700054A10Rik, AK089884, and Daw1 as the core lncRNA with high significance in diabetic cardiomyopathy.

Left ventricular hemodynamics of control and db/db mice. Table 1 lists left ventricular (LV) hemodynamic parameters of C57BL/6J and db/db mice at 6 and 20 weeks of age. There were no significant differences in heart rate, LV end-systolic pressure, LV end-diastolic pressure, dP/dt max (maximum rate of increase of LV developed pressure), and dP/dt min (maximum rate of decrease of LV developed pressure) between db/db and C57BL/6J mice at both 6 and 20 weeks of age (P > 0.05, n = 10 mice/group).
Db/db mice at 20 weeks of age developed early DCM. The dimensions and function of the LV measured by echocardiography are shown in Table 2. LV anterior and posterior walls were increased in db/db mice at both 6 and 20 weeks of age compared with age-mated C57BL/6J mice (n = 12-13 mice/group, P < 0.05). There were no significant differences in LV internal diameters at both end diastole and end systole and fractional shortening between db/db mice and age-matched C57BL/6J mice (P > 0.05). Mitral E/A ratio was comparable between db/db and C57BL/6J mice at 6 weeks of age (P > 0.05, n = 12-13 mice/group). However, it was significantly lower in 20-week-old db/db mice than C57BL/6J controls (P < 0.05). Kruskal-Wallis test followed by Dunn's test was used to analyze multiple group comparisons. *P < 0.05 versus C57BL/6J group at 6 weeks old; # P < 0.05 versus C57BL/6J group at 20 weeks old (n = 12-13 mice/group). 6   www.nature.com/scientificreports www.nature.com/scientificreports/ and 197 and 291 mRNAs were down-regulated in 6 and 20-week-old db/db mice, respectively (Fig. 6A,B) (fold change ≥ 2.0 and P < 0.05). Hierarchical clustering analysis revealed different de-regulated mRNA signature in diabetic young and older mice compared with controls ( Fig. 7A,B). Deregulated mRNAs distributed in 20 chromosomes in both 6 and 20 weeks-old db/db mice. Among 20 chromosomes, chromosome 7 had the maximum number of aberrantly expressed mRNAs (Fig. 7C,D). Compared with controls, 32 and 42 mRNAs in Bioinformatics analysis of differentially expressed mRNAs in db/db mouse hearts with and without DCM. Figures 8 and 9 show GO analysis of 120 differentially expressed mRNAs in db/db mouse hearts with and without DCM. In 6-week-old db/db mouse hearts, top 30 up-regulated mRNAs were linked with 426 GO terms in the biological process network, 28 in cellular component networks, and 90 in molecular function network (Fig. 8A). In 20-week-old mouse hearts, top 30 up-regulated were associated with 529 GO terms in biological process network, 97 in cellular component network, and 112 in molecular function network (Fig. 8B). In 6-week-old db/db mouse hearts, the down-regulated mRNAs were associated with 926 GO terms in biological processes network, 67 in cellular component network, and 91 in molecular function network (Fig. 9A). In 20-week-old db/db hearts, the down-regulated mRNAs were linked 442 GO terms in biological processes network, 36 in cellular component network, and 96 in molecular function network (Fig. 9B). Figure 10 shows the KEGG pathway analysis of deregulated mRNAs in db/db mouse hearts. In 6-week-old db/db diabetic hearts, up-regulated mRNAs were significantly enriched for renin-angiotensinogen system,  www.nature.com/scientificreports www.nature.com/scientificreports/ endocrine and other factor-regulated Ca 2+ reabsorption, and taurine and hypotaurine metabolism (Fig. 10A). In 20-week-old db/db diabetic hearts, up-regulated mRNAs were enriched for renin-angiotensinogen system, endocrine and other factor-regulated Ca 2+ reabsorption, and chemokine signaling pathway (Fig. 10B). In 6-week-old db/db diabetic hearts, down-regulated mRNAs were significantly enriched for mitogen-activated protein kinase signaling pathway, pancreatic cancer, and breast cancer (Fig. 10C). In 20-week-old db/db diabetic hearts, down-regulated mRNAs were enriched for influenza A, toxoplasmosis, and tuberculosis (Fig. 10D). www.nature.com/scientificreports www.nature.com/scientificreports/ LncRNA-mRNA co-expression networks. Significantly co-expressed lncRNAs-mRNAs (Pearson correlation > 0.995 or < −0.995 and P < 0.05) were identified and assembled into co-expression networks. There were 10,100 connections between 654 lncRNAs and 439 mRNAs in db/db mouse hearts at 6 weeks of age and 21,463 connections between 1601 lncRNAs and 1277 mRNAs in db/db mouse hearts with DCM. Figure 11 shows the co-expression network of 5 lncRNAs having the maximum number of connections with mRNAs. In 6-week-old db/db mouse hearts, top 5 lncRNAs (Gm20717, Gm14492, Gm11209, Epn2, Mapk14) were connected with 76 mRNAs (Fig. 11A). These 5 lncRNAs were correlated with the development of myocardial tissues (GO: 0055024), positive regulation of myocardial development (GO: 1901863), positive regulation of striated muscle tissue development (GO: 0045844), positive regulation of muscle organ development (GO: 0048636), tissue regeneration (GO: 0042246), and hydrolase activity (GO: 0016811) (Pearson correlation > 0.995 or < −0.995 and P < 0.005) (Fig. 12A). In 20-week-old db/db mouse hearts, top 5 network nodes (BC038927, G730013B05Rik, 2700054A10Rik, AK089884, Daw1) were linked with 126 mRNAs (Fig. 11B). These lncRNAs were correlated with action potential (GO: 0086001), membrane depolarization (GO: 0070252), conduction (GO: 0061337), contraction (GO: 0086003), and actin filament based movement (GO: 0030048) of cardiac muscles (Pearson correlation > 0.995 or < −0.995 and P < 0.005) (Fig. 12B).
Validation of lncRNAs and mRNAs using quantitative reverse transcriptional-polymerase chain reaction. To validate the RNA-sequence data, we used quantitative reverse transcriptional-polymerase chain reaction (qRT-PCR) to analyze 12 lncRNAs and 12 mRNAs in 6 and 20-week-old mouse hearts. qRT-PCR analysis www.nature.com/scientificreports www.nature.com/scientificreports/ revealed that the expression levels of the lncRNAs, AA465934, Ccdc92, AK030243, AA465934, Sart3, and Meg3, were up-regulated in db/db mouse hearts, while those of XLOC_010800, 3300005D01Rik, Gm12913, AK143388, Nxn, and AK038009 were down-regulated in db/db mouse hearts. Similarly, the expression levels of the mRNAs, Klk1b16, AI593442, Fbxw19, AI593442, Klk1b16, and Cmtr1 were up-regulated in db/db mouse hearts, whereas those of Nkd1, Ttc39a, Cd74, Gdf15, 5031410I06Rik, and Gm10220 were down-regulated in db/db mouse hearts. The qRT-PCR results of both lncRNAs and mRNAs were mostly consistent with the microarray data.

Discussion
The results of the present study demonstrate that db/db mice have diabetes, obesity, and LV hypertrophy at 6 weeks of age and develop early DCM at 20 weeks of age. The profiling of both lncRNA and mRNA expression is different between db/db mouse hearts without DCM and with early DCM. We identified 754 lncRNAs that were aberrantly expressed in db/db mouse hearts with early DCM. Among them, the lncRNAs, BC038927, G730013B05Rik, 2700054A10Rik, AK089884, and Daw1, have maximum connections with differentially expressed mRNAs. Bioinformatics analysis reveals that these 5 lncRNAs are closely associated with membrane depolarization, action potential conduction, contraction, and actin filament-based movement of cardiac cells. These results suggest that BC038927, G730013B05Rik, 2700054A10Rik, AK089884, and Daw1, may be the core lncRNA with high significance.
The db/db mouse is the most widely used rodent models in obesity-induced T2DM research 18,19 . Db/db mice suffer from a deficiency or mutation of leptin receptors on chromosome 4 20 . Leptin is an adipose tissue hormone and functions to regulate food intake, energy homeostasis, thermogenesis, reproduction, hematopoiesis, skeletal growth, and neuroprotection by activating the cognate leptin receptors 21 . In the db/db mouse, leptin fails to exert its actions due to the deficiency or mutation of leptin receptors, displaying hyperglycemia, obesity, hyperinsulinemia, and dyslipidemia 17 . This was confirmed by our present study showing that the db/db mice at both 6 and 20 weeks of age had increased body weight, blood glucose, and plasma insulin levels.
In humans and monkeys, DCM is characteristic of cardiac (both diastolic and later systolic) dysfunction that occurs independently of a recognized cause such as coronary artery disease or hypertension 8,22,23 . In the present study, db/db mice at 6 weeks of age had normal cardiac function, despite hyperglycemia, obesity, and cardiac hypertrophy. These results are consistent with our previous report in db/db mice with the genetic background of C57BLKS/J (C57BL/KsJ-lepr db /lepr db ) 18 . Intriguingly, the db/db mice at 20 weeks of age had normal systolic function (fractional shortening) but decreased diastolic function (mitral E/A ratio). This diastolic dysfunction was independent of LV pressure. Moreover, the electrocardiogram examination of 20-week-old db/db mice revealed that there was no myocardial ischemia (data not shown). Thus, depressed diastolic function in 20-week-old db/ db mice can be attributed to early DCM.
A growing number of studies indicate that lncRNAs play crucial roles in various cardiovascular diseases, including cardiac hypertrophy, heart failure, and diabetic vascular complications [24][25][26][27][28] . Recent studies reported that a few lncRNAs, including H19, metastasis-associated lung adenocarcinoma transcript 1 (MALAT1), and myocardial infarction-associated transcript (MIAT), were involved in the pathogenesis of DCM in type 1 diabetic animals 16,29,30 . We identified 754 lncRNAs that were differentially expressed in diabetic db/db mouse hearts with early DCM. It is likely that multiple deregulated lncRNAs contribute to the pathogenesis of early DCM in the T2DM db/db mice.
In the present study, 754 deregulated lncRNAs in db/db mouse hearts with early DCM were distributed on 21 chromosomes. However, these lncRNAs were not equally distributed across all chromosomes. Compared with other chromosomes, chromosome 2 had the maximum number of deregulated lncRNAs in db/db mouse hearts with and without DCM. Diabetic cardiovascular disease is increasingly identified to associate with genetic susceptibility [31][32][33][34] . Therefore, chromosome 2 may be more likely to carry lncRNAs susceptible to the pathology of DCM.
LncRNAs regulate epigenetically the expression of broad ranges of genes in cardiomyocytes 13,35 . To gain a better insight into the biological function of deregulated lncRNAs and conduct the analysis of transcribed sequences, it is advantageous to arrange them into several clusters 4,36 . Based on the relationship between lncRNAs and their affiliated protein-coding genes, lncRNAs detected by Arraystar Array are characterized as antisense, bidirectional, intronic, sense overlapping and intergenic lncRNAs 36 . The majority of differentially expressed lncRNAs in our study were intergenic (50%) and antisense (30%), which contributes to more than three-quarters of the www.nature.com/scientificreports www.nature.com/scientificreports/ deregulated lncRNAs. LncRNAs are transcribed from regions of at least 5 kb, from protein-coding genes. They can modulate the expression of target genes, and the target genes can be scattered across the genome. Antisense lncR-NAs are capable of altering the expression of protein-coding genes and affect their protein-coding counterparts www.nature.com/scientificreports www.nature.com/scientificreports/ through various mechanisms 37 . The many abnormally expressed intergenic lncRNAs and antisense lncRNAs in DCM indicate that lncRNAs may regulate protein-coding genes throughout the development of DCM.
Our microarray data showed that a large number of mRNAs were aberrantly regulated in db/db mouse hearts with and without DCM. Altogether, 447 and 662 mRNAs were found to differentially express in db/db hearts without DCM and db/db hearts with DCM, respectively. Overall, more mRNAs were up-regulated than down-regulated. Deregulated mRNAs were distributed on 21 chromosomes, among which chromosome 7 had the maximum number of aberrantly expressed mRNAs.
LncRNAs are increasingly identified to play roles in a broad variety of principal biologic activities 35,38 . The function of lncRNAs proximately associates with their correlated protein-coding genes, thus studying the genomic background of lncRNAs can help predict their essential biological function 39 . Based on the GO and KEGG database analyses of lncRNAs, protein-coding genes, and mRNAs, up-regulated genes were significantly enriched for renin-angiotensinogen system, endocrine and other factor-regulated Ca 2+ reabsorption, and taurine and hypotaurine metabolism in db/db diabetic hearts without DCM and for renin-angiotensinogen system, endocrine and other factor-regulated Ca 2+ reabsorption, and chemokine signaling pathway in db/db diabetic hearts with DCM. Whereas down-regulated mRNAs were significantly enriched for mitogen-activated pancreatic cancer, protein kinase signaling pathway, and breast cancer in db/db diabetic hearts without DCM and for influenza A, toxoplasmosis, and tuberculosis in db/db diabetic hearts with DCM. Previous studies have shown that cancers, viral infections, and tuberculosis are higher in diabetic patients [40][41][42][43] . Thus, these down-regulated mRNAs in diabetes may be involved in increased cancers and infection.
The integrated analysis of the lncRNA-mRNA co-expression networks revealed that 5 deregulated lncR-NAs having maximum connections with differentially expressed genes were BC038927, G730013B05Rik, 2700054A10Rik, AK089884, and Daw1. Functional analysis identified these 5 lncRNAs were correlated with membrane depolarization of cardiac muscles like actin mediated cell contraction (GO: 0070252), cardiac muscle cell contraction (GO: 0086003), cardiac conduction (GO: 0061337), cardiac muscle cell action potential (GO: 0086001), and actin filament based movement (GO: 0030048). To our knowledge, there are no previous studies reporting the function and role of these lncRNAs in DCM. However, our results suggest that these 5 lncRNAs may be core lncRNA with high significance for early DCM.
In summary, the results of our investigation indicate that db/db mice at 20 weeks of age develop early DCM. Among aberrant expression of the 754 lncRNAs that are associated with the pathogenesis of early DCM, BC038927, G730013B05Rik, 2700054A10Rik, AK089884, and Daw1, have maximum associations with mRNAs. Given that these 5 lncRNAs are associated with membrane depolarization, action potential conduction, contraction, and actin filament-based movement of cardiac cells, they may serve as important therapeutic targets for DCM.

Materials and Methods
Animals. Obese male db/db and C57BL/6J control mice were purchased from The Jackson Laboratory (Bar Harbor, ME, USA). The animals were kept on a 12-h light-dark cycle in a temperature-controlled room. The animal care and all experimental procedures were performed in accordance with the NIH Guide for the Care and Use of Laboratory Animals (Institute for Laboratory Animal Research, National Academy of Sciences, 8th edition, 2011), and experimental protocols were approved by the Institutional Animal Care and Use Committee (IACUC) at the Medical College of Wisconsin (Milwaukee, WI, USA).
Our pilot experiments showed the db/db mouse had the normal systolic function of the left ventricle but started to have impaired diastolic function at 20 weeks of age. Based on the preliminary studies, the db/db mice were divided into 2 experimental groups: diabetic mice without DCM (db/db mice at 6 weeks of age) and diabetic mice with early DCM (db/db mice at 20 weeks of age). Age-and gender-matched C57BL/6J mice were used as controls.
Measurements of blood glucose and insulin. C57BL/6J and db/db mice at either 6 or 20 weeks of age were fasted for 6 h (12-13 mice/group). Under the anesthesia of 80 mg/kg pentobarbital, a thoracotomy was performed, and the left ventricle was punctured with a 27 gauge needle, as described previously 44 . Blood glucose was measured with a blood gas analyzer (ABL-725 Radiometer, Radiometer America Inc., Westlake, OH, USA) 45 . In ice-chilled heparinized tubes, plasma was immediately separated and stored frozen at −80 °C. Plasma insulin levels were measured by radioimmunoassay in 20 μl aliquots of plasma using a commercial kit (Linco Research Inc., St Charles, MI, USA) 46 .
LV hemodynamic measurement. C57BL/6J and db/db mice were anesthetized by the inhalation of 2.0% isoflurane and oxygen (n = 10 mice/group). The mice were ventilated with room air supplemented with 100% oxygen at approximately 102 breaths/min, as described 44,47 . The right carotid artery was cannulated with a Millar Mikro-Tip Pressure Transducer Catheter (1.4-Fr, model SPR 671, Millar Instruments, Inc.; Houston, TX, USA), and the catheter was placed in the middle of the LV chamber to measure left ventricular systolic and diastolic pressure, as described 48 . Before every measurement of left ventricular pressure (LVP), the catheter was calibrated electronically in vitro by submerging its tip in the saline solution of 37 °C and applying an external pressure of 100 mmHg. The catheter was connected to an ADInstrument pressure transducer (MLT0380/D, ADInstruments, Colorado Springs, CO, USA) and a Powerlab data acquisition system (ADInstruments). After a 30 min of stabilization, blood pressure was continuously recorded for 20 min 48 . The LVP signal was monitored, and dP/dt max and dP/dt min were determined. Body temperature was maintained between 36.8 °C and 37.3 °C throughout the experiment by using a heating pad (Model TC-1000, CWE Inc.; Ardmore, PA, USA). (2019) 9:15345 | https://doi.org/10.1038/s41598-019-51872-9 www.nature.com/scientificreports www.nature.com/scientificreports/ Transthoracic echocardiography. C57BL/6J mice and db/db mice were sedated by the inhalation of 1.50% isoflurane and oxygen (n = 12-13 mice/group). Non-invasive transthoracic echocardiography was performed with a VisualSonics Vevo 770 High-resolution Imaging System (Toronto, Canada) equipped with a 30 MHz transducer (Scanhead RMV 707), as described 18,49 . LV dimensions and fractional shortening were measured by two-dimension guided M-mode method. Pulsed Doppler waveforms recorded in the apical-4-chamber view were used for the measurements of the peak velocities of mitral E (early mitral inflow) and A (late mitral inflow) waves.

RNA extraction.
After echocardiographic examination was completed, the C57BL/6J and db/db mice were euthanized to harvest the LV (n = 12-13 mice/group). Myocardial tissues were immersed immediately into liquid nitrogen, and total RNAs were isolated using Triazol Reagent (Invitrogen, Carlsbad, CA, USA), as described 50,51 . The extracts were further quantified using Nano Drop ND-1000. Presence of DNA contaminants was assessed by agarose gel electrophoresis on a denaturing gel, whereas RNA integrity was determined using Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Array hybridization was performed, according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technologies) with minor modifications. Briefly, mRNA was purified from total RNAs after removal of rRNA with mRNA-ONLY ™ Eukaryotic mRNA Isolation Kit (Epicentre Biotechnologies, Madison, WI, USA). Each sample was amplified and transcribed into fluorescent cRNA along the entire length of the transcripts without 3ʹ bias utilizing a mixture of oligo(dT) and random priming method with Arraystar Flash RNA Labeling Kit (Arraystar Inc., Rockville, MD, USA). The labeled cRNA was purified by RNeasy Mini Kit (Qiagen, Valencia, CA, USA). The concentrations and specific activity of the labeled cRNA (pmol Cy3/μg cRNA) were measured by Nano Drop ND-1000. One μg of each labeled cRNA was fragmented by adding 5 μl 10× blocking agent and 1 μl of 25× Fragmentation Buffer, then heated the mixture at 60 °C for 30 min, finally 25 μl 2 × GE Hybridization buffer was added to dilute the labeled cRNA. Fifty μl of hybridization solution was dispensed into the gasket slide and assembled to the lncRNA expression microarray slide. The slides were incubated for 17 h at 65 °C in an Agilent hybridization oven. The hybridized arrays were washed, fixed, and scanned by the Agilent DNA Microarray Scanner (Agilent Technologies).
Microarray and bioinformatics analysis. The profiling of lncRNA and mRNA expression was performed using Arraystar Mouse LncRNA Microarray V3.0 (Arraystar Inc.), which detected 35,923 lncRNAs and 24,881 mRNAs 52 . Agilent Feature Extraction software version 11.0.1.1 (Agilent Technologies) was used to analyze acquired array images. Quantile normalization and subsequent data processing were performed using the GeneSpring GX v12.1 software package (Agilent Technologies). After quantile normalization of the raw data, lncRNAs and mRNAs that at least 4 out of 16 samples had flags in Present or Marginal ("All Targets Value") were chosen for further data analysis. Statistically significant changes between the two groups were considered when the fold changes for differentially expressed lncRNAs and mRNAs were larger than 2.0, and the P value for t-test was less than 0.05. Differentially expressed lncRNAs and mRNAs with statistical significance were identified through Volcano Plot filtering between two groups. Hierarchical clustering was performed using the R software (version 2.15). GO (www. geneontology.org) and KEGG enrichment analyses were made to annotate the potential functions of differentially expressed lncRNAs. To identify the potential transcriptional regulators, the analysis of upstream regulators was conducted using the Ingenuity Pathway Analysis software (Ingenuity Systems Inc., Redwood City, CA, USA). Only transcriptional regulators that appeared up-or down-regulated in the microarray study were considered for the analysis.
LncRNA-mRNA correlation analysis. To reveal potential association of the differentially expressed lncRNAs with mRNAs in DCM, we built the lncRNA-mRNA co-expression networks using Cytoscape software (v3.4.0), according to the normalized signal intensity of individual genes 53 . Briefly, microarray data were pre-processed by using the average expression value of all transcripts expressed from the same gene (both mRNA and lncRNA). The data were then screened for differentially expressed lncRNAs and mRNAs whose expression levels positively or negatively correlated. For each pair of lncRNA-mRNA, Pearson correlation test was conducted to detect significant correlation. Only strongly correlated (r 2 ≥ −0.9, P < −0.01) pairs were used to construct the networks and generate visual representations. In these representations, each gene corresponded to a node, and the color of nodes represented the up-regulated or down-regulated expression of the specific gene in the microarray data.
Total RNA from the LV was extracted using Triazol Reagent, as described above. Chloroform was added, and samples were centrifuged to facilitate phase separation. The aqueous phase was extracted and combined with ethanol in miRNeasy Mini spin columns (Qiagen). Total RNA was eluted in RNase-free water, and the concentration of extracted total RNA was quantified by the Epoch spectrophotometer (Biotek, Winooski, VT). Samples were considered pure if the A260/280 ratio was between 1.9 and 2.0. One µg of total RNA from each sample was used to generate cDNA using miScript Reverse transcriptase mix, nucleic mix, and HiFlex Buffer (Qiagen). To analyze the lncRNA expression, a master mix (25 μl/well) containing the template cDNA (4.5 ng/well), RNase-free