Whole transcriptome RNA-seq reveals key regulatory factors involved in type 2 diabetes pathology in peripheral fat of Asian Indians

The prevalence of Type 2 Diabetes has reached an epidemic proportion particularly in south Asian countries. We have earlier shown that the anatomical fat distribution, termed ‘thin fat phenotype’ in this population indeed plays a major role for their T2D-predisposition it is indeed the sick fat or adiposopathy, which is the root cause of metabolic syndrome and diabetes and affects both—peripheral, as well as visceral adipose tissue compartments. In present study, we have attempted to unravel the altered regulatory mechanisms at the level of transcription factors, and miRNAs those may likely accounts to T2D pathophysiology in femoral subcutaneous adipose tissue. We prioritized transcription factors and protein kinases as likely upstream regulators of obtained differentially expressed genes in this RNA-seq study. An inferred network of these upstream regulators was then derived and the role of TFs and miRNAs in T2D pathophysiology was explored. In conclusions, this RNS-Seq study finds that peripheral subcutaneous adipose tissue among Asian Indians show pathology characterized by altered lipid, glucose and protein metabolism, adipogenesis defect and inflammation. A network of regulatory transcription factors, protein kinases and microRNAs have been imputed which converge on the process of adipogenesis. As the majority of these genes also showed altered expression in diabetics and some of them are also circulatory, therefore they deserve further investigation for potential clinical diagnostic and therapeutic applications.

www.nature.com/scientificreports/ etc. These results strengthen the belief that it is indeed the sick fat or adiposopathy, which is the root cause of metabolic syndrome and diabetes and affects both-peripheral, as well as visceral adipose tissue compartments 6 . We further attempted to link the molecular basis of aforementioned thin fat phenotype in Asian Indian diabetics with that of a monogenetic disorder, lipodystrophy which shows marked similarity with T2D in clinical presentation 7 . A significant overlaps between physical and functional protein-protein interaction networks of differentially expressed genes (DEGs) in peripheral fat of diabetics and annotated lipodystrophy genes was found. To some extent, we could be able to stratify patients in terms of differentially expressed lipodytrophy genes and henceforth opined that functionally disturbed expression of lipodystrophy genes might play a role in T2D pathogenesis.
Transcription factors (TFs) are proteins that bind to the DNA and help initiate a program of increased or decreased gene transcription. Many human diseases have been associated with mutations in transcription factors for example mutation in hepatocyte nuclear factors (HNFs) results in a rare form of diabetes called MODY (Maturity onset diabetes of the young) 8 . As therapeutic modulation of TFs by drugs may lead to dramatic changes in overall gene expression pattern and could ameliorate the disease condition, identification of key TFs in T2D pathophysiology therefore promise toward the development of better anti-diabetic agents.
microRNAs are short noncoding RNAs which degrade those mRNAs whose product are not required in normal cellular physiology. Conversely, an altered miRNA expression pattern leads to a derangement in cellular physiology. Various human microRNA-disease associations have been validated experimentally and this knowledge can help us in the functional interpretation of RNA-seq datasets.
In present study, we have attempted to unravel the altered regulatory mechanisms at the level of transcription factors, and miRNAs those may likely accounts to T2D pathophysiology in femoral subcutaneous adipose tissue. We conducted whole transcriptome RNA-seq on a subset of non-diabetic, and diabetic peripheral subcutaneous adipose tissue samples which we have earlier used in our gene expression microarray study 6 . We have also prioritized transcription factors and protein kinases as likely upstream regulators of obtained differentially expressed genes in RNA-seq study. An inferred network of these upstream regulators was derived and the role of TFs and miRNAs in T2D pathophysiology was explored in the light of existing literature.
Protein kinases (PKs) have become the second most important group of drug targets after G-protein coupled receptors as deregulation in majority of signal transduction pathways have been found to be associated with abnormal phophorylation by kinases 9 . We therefore attempted to derive a minimal molecular pathway from the inferred network, comprising top-ranked TFs, PKs, lipodystrophy genes, and the central mediator of adipogenesis-PPARγ. Three kinases were identified using network analysis and proposed as potential drug targets.

Results
A total of 2752 differentially expressed (DE) genes were obtained by DESeq2 10 (P < 0.01; logFC ± 2) between the femoral fat of non-diabetics and diabetics (See Supplementary File S1). Figure 1 shows qPCR results for the ten selected genes that validated the directionality of gene expression obtained from RNA-seq (See Supplementary File S2).  www.nature.com/scientificreports/ Top 2000 genes were then selected on the basis of their p-value for subsequent pathway enrichment analysis against the human collection of KEGG pathways 11 using STRING db 12 as it supports only a network less than 2000 nodes. We found pathways linked with glucose-and fat-metabolism (shown in bold text in Table 1) which were in line with our earlier enrichment results for microarray datasets of femoral fat.
We also found a total of 65 miRNAs to be differentially expressed between non-diabetic and diabetic groups (P < 0.05). Eight of these miRNAs have been found to be associated with type 2 diabetes, diabetic retinopathy, nephropathy, obesity, and inflammation (Table 2).
CyTargetLinker 13 analysis further revealed that these miRNAs were regulating expression of 214 differentially expressed genes (Fig. 2). mir133b has been reported to attenuate GLUT4-mediated glucose update by regulating expression of its regulator Kruppel-like transcription factor 15 (KLF15) 14 . mir210 also attenuates insulin signalling by targeting PTPN1 and it's up-regulation has been reported in T1D patients 15 . miR-22-3p activates gluconeogenesis by down-regulating the expression of transcription factor 7 (TCF7) which is a negative regulator of enzymes of the gluconeogenic pathway 16 . Down-regulation of miR-25 has also been observed across multiple studies in diabetic nephropathy 17 . Up-regulation of miR-27a has been reported in rat model of T2D 18 as well as a its high circulating level has been reported in patients with metabolic syndrome and T2DM 19 . miR-27a has been observed a negative regulator of peroxisome proliferator-activated receptor-γ (PPAR-γ) and CCAAT/enhancer-binding protein-α, which are extremely important regulators of adipogenesis 20 . Clinically its expression was found to be positively correlated with fasting glucose level 19 . The involvement of miR-27b in diabetes nephropathy has been reported. It targets pathways related to ER stress and chronic ER stress has previously been identified in peripheral blood mononuclear cells (PBMCs) from T2DM patients 21 . Lin et al. 20 has found a reduced expression of miR-27a www.nature.com/scientificreports/ and -b during adipogenic differentiation of 3T3-L1 cells in a microarray study, which was correlated with an increase of expression of PPARγ. miRNA-27b also contributes to lipopolysaccharide-mediated PPARγ mRNA destabilisation 22 . miR-30 has also been found down-regulated in patients with type 2 diabetes 23 . Deregulation of miRNA-503 has also been reported to contribute to T2D-induced impairment of endothelial function 24 . It also shows differential expression between T2D from obese-T2D patients 25 . All the diabetes-specific miRNAs were circulatory and therefore could be further followed up for their potential use as biomarkers. We also found differential expression of other miRNAs such as hsa-mir-6784, hsa-mir-4707, hsa-mir-647, hsa-mir-3652, hsa-mir-6757, hsa-mir-1-1, hsa-mir-3671, hsa-mir-6505, hsa-mir-1229, and hsa-mir-LeT7F1. The deregulated expression of hsa-mir-1-1 has been found in myocardial infarction and other cardiac complications. Similarly, hsa-mir-647 was also associated with ischemic stroke. The expression of hsa-mir-1229 was linked with Colorectal Carcinoma. The down regulation of hsa-mir-LeT7F1 has been reported in many cancers. eXpression2Kinases (X2K) 26 constructed a network of 110 proteins (Fig. 3) for the inputted list of differentially expressed genes. This network comprised the likely upstream TFs and PKs regulating DE genes.
Top 10 transcription factors in the network are presented in Table 3. All these transcription factors were found to be involved in processes associated with adipogenesis, or glucose homeostasis and seven of them were also differentially expressed (P < 0.05). Polycomb Repressive Complex 2 (PRC2) is involved in epigenetic silencing of PPARγ and two enriched transcription factors SUZ12 and EZH2, as part of PRC2 negatively regulate its transcriptional activation 27 . Similarly, Hepatocyte nuclear factor 4α (HNF4α) together with its target Hes6 maintains low expression level of PPARγ in liver 28 so its expression in peripheral fat reflect poor adipogenesis. MyoD up-regulates PPARγ expression and has been reported to promote adipocyte differentiation 29 . ESR1 (Estrogen Receptor 1) has been reported to physically associates with PPARγ and functionally interferes with its signalling 30 . The roles of Nuclear factor erythroid 2-related factor 2 (NRF2) in the regulation of metabolism, inflammation, mitochondrial physiology, and immune responses has been recognised 31 and its modulation is extensively being studied for treatment of diseases that are caused by oxidative stress 32 . NRF2-knockout mice were found to be partially protected from high fat diet-induced obesity and developed a less insulin-resistant phenotype 33 . Transcription factor EGR1 has also been reported to act as a negative regulator of 3T3-L1 adipocyte differentiation 34 . FOXA2 is a negative regulator of adipocyte differentiation. It enhances the insulin receptor (IR), insulin receptor substrate-2 (IRS2), glucose transporter-4 (GLUT4), HK2, and UCP2 and UCP-3 genes in preadipocytes and adipocytes 35 . Another transcription factor REST has been reported to negatively regulate the expression of several neuronal genes, such as glutamate receptor 7 and synuclein γ and expression of these genes in adipocytes has been reported which suggests the role of REST in adipocyte differentiation and function in both white as well as brown adipose tissues 36 . The role of TP65 is adipose biology is known and PPARγ agonists, including thiazolidinediones (TZDs) have been reported to induce anti-proliferation, differentiation, and apoptosis in adipocytes through inducing expression of endogenous or exogenous TP63 37 .
Three lipodystrophy-associated genes, Histone deacetylase 3 (HDAC3), Nuclear receptor subfamily 0 group B member 2 (NR0B2), and small ubiquitin-like modifier 1 (SUMO1) were also present in this network. HDAC3 induces lipodystrophy by inhibiting glyceroneogenesis in adipocytes through repressing cytosolic phosphoenolpyruvate carboxykinase (PEPCK) 38 . NR0B2 is transcriptionally regulated by PPARγ and has been found to be associated with insulin resistance 39 . SUMO1 has been found to repress transcriptional activity of PPARγ through its sumoylation and consequently impairing its pro-adipogenic function 40 .
The gene HNF4A (Hepatocyte Nuclear Factor 4 Alpha) has been associated with monogenic autosomal dominant non-insulin-dependent diabetes mellitus type I and has been identified in GWAS studies in people of south Indian ancestrary alongwith another gene TMEM163 (Transmembrane Protein 163) 41 . Its variants also showed association with decreased fasting plasma insulin and homeostatic model assessment of insulin resistance, indicating a plausible effect through impaired insulin secretion 42 . www.nature.com/scientificreports/ A dual metric approach was then used to calculate an average score of node using degree-score, and bottleneck-score that reflects the importance of each node in this minimal network. Top 7 nodes are presented in Table 5. Three kinases: MAPK3, MAPK14, and AKT1 were amongst these nodes.
All these kinases have been found to be associated with diabetes-related KEGG pathways-AGE-RAGE signaling pathway in diabetic complications (hsa04933), and Endocrine resistance (hsa01522). MAPK3, and AKT1 were involved in Insulin signaling pathway (hsa04910), and Type II diabetes mellitus (hsa04930). AKT1 is also involved in Adipocytokine signaling pathway (hsa04920), Regulation of lipolysis in adipocytes (hsa04923), Insulin resistance (hsa04931) pathways. Clearly their deregulation seems to play role in the T2D pathogenesis.
AKT Inhibitors are considered as very promising drug target in cancer research as dysregulation in phosphoinositide 3-kinase/AKT/mammalian target of rapamycin pathways have been found to associate with many     www.nature.com/scientificreports/ cancers. However, due to connection of these pathways with insulin signalling, AKT Inhibitors have been reported to induce hyperglycaemia 43 so their use as drug target in diabetic complications seems not appropriate. Inhibition of MAPK3 (ERK1) inhibitor in hypertrophic 3T3-L1 adipocytes has been reported to ameliorate the dysregulation of adipocytokine expression and suppressed the enhanced lipolysis activity. It has also shown to control blood glucose level in the diabetic mice 44 . A genome-wide association study suggests that MAPK14 is associated with diabetic foot ulcers 45 . MAPK14 (p38) is a stress-activated protein kinases (SAPK) and is involved in biological functions, such as inflammation, differentiation, proliferation, and apoptosis 46 . Pentosan polysulfate-an inhibitor of p38 MAPK pathway has been found to attenuate apoptosis and inflammation by down-regulating this pathway in high glucose treated human renal proximal tubular epithelial cells 47 .

Discussion
In summary, this study find that in diabetics peripheral subcutaneous adipose tissue show molecular pathology characterized by differential expression of genes enriching glucose, lipid and protein metabolic pathways as well as adipogenesis, inflammation and coagulation related pathways. The imputed regulatory network of this transcriptome contained several transcription factors and protein kinases which converge on the process of adipogenesis. The majority of them were also differentially expressed in diabetics. Out of 52 differentially expressed miRNAs, eight are circulatory and are known to be associated with diabetes and its complications. Additionally few like miR-27a are also known to regulate the process of adipogenesis. Interestingly many of the regulatory factories were mapped to genetic loci previously showed association with diabetes and related intermediate phenotypic traits in genome wide association studies (GWAS) In other words these finding shed light on the regulatory mechanism and genome to phenome pathways of diabetes in peripheral subcutaneous adipose tissue. Peripheral adipose tissue is generally considered protective in the pathogenesis of insulin resistance, but not so in the case of Asian Indians diabetics, who show thin fat phenotype. Recently we have reported that in diabetics peripheral subcutaneous adipose tissue on transcription profiling, show not only molecular pathology consistent with the so called "sick fat" or "adiposopathy", but also its modules of co-expressed genes showed an association with several disease related intermediate phenotypic traits. In other words, these findings suggested the role of this adipose tissue pathology in the pathogenesis of insulin resistance and diabetes. Interestingly on functional analysis these modules of co-expressed genes enriched several adipogenesis and inflammation related pathways. These findings are consistent with the theory that it is the limitation of adipogenesis in subcutaneous adipose tissue and inability to expand and store excess fat lead to ectopic fat deposition in liver and muscle. The ectopic fat deposition in liver and muscle play an important role in the pathogenesis of insulin resistance and diabetes. The results of present study add to the current knowledge in this field by identification of network of transcription factors, protein kinases and miRNAs regulating this molecular pathology.
It was observed in this study that eight microRNAs known to be associated with T2D and its co-morbidities like diabetes retinopathy, and nephropathy were differentially expressed (P < 0.05) in diabetics.Over 214 DE genes, accounting to ~ 7% of all DE genes were being regulated by these microRNAs at stringent P < 0.01. Several of them in previous studies were found to be associated with insulin signalling, inflammation and adipogenesis, thereby suggesting their role, not only in sick fat or adiposopathy, but possibly also in the aetiology of cardiovascular and microvascular complications of diabetes. As most of them are circulatory miRNAs, they can not only serve as targets of new drug discovery, but also the biomarkers of adiposopathy. As a clinical indicator of adiposopathy such biomarkers can be better defined than the presently employed indicator-metabolic syndrome. Therefore, they deserve further investigations towards potential clinical diagnostic and therapeutic application.
Another important finding of this study is that a small network of 110 proteins, comprising upstream regulators of differentially expressed genes (i.e. transcription factors, protein kinases and intermediate proteins connecting them) was imputed on the basis of the differentially expressed genes. Interestingly more them 50% of them were also differentially expressed (P < 0.05). These results affirm the marked biological significance of this RNA-seq expression profiles. Interestingly, almost all of the top enriched upstream transcription factors are known to regulate genes implicated in the process of adipogenesis. Moreover, several of these genes had also shown linkage disequilibrium with diabetes related traits in previously published genome wide association studies. Therefore, these findings suggests that some of diabetes susceptibility genes possibly manifest clinically by influencing this regulatory network and control the process of adipogenesis in peripheral subcutaneous adipose tissue. Adipogenesis defects in peripheral subcutaneous adipose tissue can explain various aspects of "thin fat" phenotype like relatively thiner limbs (i.e. smaller mid arm circumference) in Asian Indian newborns, lesser fat mass in limbs of diabetics as compared to non-diabetics and relatively high insulin resistance for given body mass index in them etc 47 . Interestingly we recently observed that insulin resistance in diabetics showed positive correlation with adipocyte size (a parameter of poor adipogenesis) in the peripheral subcutaneous adipose depot, but not in the abdominal depots, thereby suggesting the important role of poor adipogenesis in the former adipose depot in the pathophysiology of diabetes in this population 49 . Identification of regulatory transcription factors, protein kinases and miRNAs in this study opens up the avenue for identification of thin fat phenotype specific biomarkers and molecular targets for drug development. In this regulatory network three lipodystrophy-associated genes were also present. As lipodystrophy is a monogenic disorder with Mendelian mode of transmission, whereas T2D is a complex multifactorial disease and its clinical picture is created by the interaction of several environmental and genetic factors such as frequent polymorphisms of multiple genes, still both exhibit similar pathophysiological alteration such as heightened insulin resistance, loss of peripheral fat, ectopic fat deposition, and adipocytokinemia etc. Presence of these genes therefore points toward the fact that thin fat phenotype might be a functional variant of lipodystrophy.
There are several potential clinical implications of the findings of this study. The higher diabetes risk among Asian Indians is generally attributed to central obesity, whereas this study find pathologic changes, particularly www.nature.com/scientificreports/ defects in adipogenesis in femoral adipose tissue in diabetics. Therefore, thigh related anthropometric parameter deserve further investigations as surrogate markers of adipose tissue pathology underlying diabetes mellitus in population. Moreover, targeting adipogenesis defects as a therapeutic strategy for the management of diabetes in population with drugs like pioglitazone is supported by the results of this study. Delineation of peripheral adipose tissue genome to phenome mechanism and pathways of insulin resistance support the concept of "adiposopathy" or "sick fat" as a formal disease for Asian Indians thin fat phenotype. The adiposopathy model of disease in addition to current gluco-centric diagnostic system has potential advantage of providing diagnostic frame for primordial prevention of diabetes.
There are several limitations of this study like relatively smaller sample size and only transcriptomic profile was generated. In order to precisely enlist potential biomarkers and novel drug targets, there is need of multi-omic investigations deciphering adipose tissue specific genome to phenome pathways and mechanisms of diabetes and insulin resistance in individuals showing thin fat phenotype. However, findings of this study, which to the best of our knowledge is the first report on regulatory mechanisms of peripheral adipose tissue dysfunction among Asian Indians, shed light on the pathways of genome to phenome correlation in diabetics.
In conclusions, this RNS-Seq study finds that peripheral subcutaneous adipose tissue among Asian Indians show pathology characterized by altered lipid, glucose and protein metabolism, adipogenesis defect and inflammation. A network of regulatory transcription factors, protein kinases and microRNAs have been imputed which converge on the process of adipogenesis. As the majority of these genes also showed altered expression in diabetics and some of them are also circulatory, therefore they deserve further investigation for potential clinical diagnostic and therapeutic applications.

Methods
The study was conducted at Sawai Man Singh Medical College, Jaipur, India. Gluteofemoral fat biopsies were obtained from five normal glucose tolerants (NGT) and five T2D individuals. The study was approved by the ethical committee of SMS Medical College, Jaipur, India, and funded by the Indian Council of Medical Research (ICMR), New Delhi. Written informed consent was obtained from the participants prior to study. All methods were performed in accordance with the relevant guidelines and regulations. The inclusion and exclusion criteria were defined in Table 6.
Total RNA was isolated from tissue biopsies using Qiagen RNeasy mini kit using manufacturer's protocol. This procedure also captured few miRNA, tRNAs, rRNAs, lncRNAs in addition to mRNAs. QC of the extracted total RNA samples was carried out on Agilent Bioanalyzer and nanodrop. mRNA paired end libraries were generated from each of the RNA samples using Illumina TruSeq mRNA library preparation kit. Sequencing was carried out on Illumina HiSeq 2500 using paired end chemistry of 2x100 bp read length. Total RNA samples with RIN value >7 were used for paired end library preparation using the TruSeq™ RNA Sample Preparation Kit version 3.0 (Illumina, Inc.) according to the manufacturer's instructions. Briefly, the Poly-A containing mRNA was purified from 2 μg of total RNA using oligo (dT) magnetic beads and was fragmented into 200-500 bp pieces in presence of divalent cations at 94°C for 5 min using an ultrasonicator. The cleaved RNA fragments were copied into first-strand cDNA using SuperScript-II reverse transcriptase (Life Technologies, Inc.) and random primers. After second-strand cDNA synthesis, fragments were end-repaired, A-tailed and indexed adapters were ligated. The products were purified and enriched with PCR to create the final cDNA library. The tagged cDNA libraries were pooled in equal ratios and used for 2Å ~ 100 bp paired-ends sequencing on a single lane of the Illumina HiSeq2500. After sequencing, the samples were de-multiplexed and indexed-adapter sequences were trimmed and other low quality bases were filtered or trimmed using in-house Perl scripts. Thus filtered, high quality reads (20x) were used for further analysis. The QC of the raw data was done using FastQC tool (http:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/). A total of 4-8.5 GB data was generated for each of the samples.
The reads were then mapped to human reference sequence genome GRCh38.84 build using HISAT2 50 with default parameters. Mapped reads were assembled into transcripts and quantified using StringTie 51 and the annotated gtf file of GRCh38.84 build. Further all the mapped & quantified results of each samples were merged and abundances of each individual genes and transcripts were calculated using StringTie. Differential expression was calculated using DESeq2 52 algorithms.
To check the reliability of RNA-seq analyses, qPCR (Real Time PCR) was done on ten genes-AKT2, BCL6, LIPE, METTL9, CD59, LMNA, FOS, CCL2, MUSTN1, and ZNF638 in thigh (femoral) subcutaneous adipose tissue of four diabetic and four non-diabetic controls subjects in triplicates. The genes were randomly selected based on their important role in diabetes, insulin resistance, lipodystrophy, lipid storage and metabolism and in inflammatory pathways. The details are attached in supplementary sheet (S2). The oligonucleotide primer sequences used were raised using freely available online tools and synthesized as per appropriate melting temperature (Tm), www.nature.com/scientificreports/ GC%. A total of 2 μg of RNA was isolated from the thigh (femoral) subcutaneous adipose tissue biopsy using Qiagen RNeasy Mini Kit (Cat No. 74104) and quantified using the microfluidic-based capillary electrophoresis system (Bio-Rad Experion). cDNA Synthesis was done using QuantiNova Reverse Transcription Kit (Cat No. 205411) and real time PCR was done using Quanti Nova SYBR Green PCR Kit (Cat. No 208052). For identification of enriched KEGG pathways by differentially expressed genes, STRING database (version 11.0) was used. To identify upstream transcription factors those regulating these differentially expressed genes, we used eXpression2Kinases (X2K) which derived an inferred network of transcription factors, protein kinases and intermediate proteins connecting them. A minimal-molecular pathway was derived from this network on the basis of connectivity of top-ranked TFs, PKs, PPARγ, and lipodystrophy genes and analyzed for therapeutic cues using network-based centrality measures like degree, and bottleneck using Cytohubba 53 .
For identification of differentially expressed miRNAs associated with diabetes, and related co-morbidities such as retinopathy, nephropathy, obesity, and inflammation, we searched Human MicroRNA Disease Database (HMDD v3.0) 54 which manually collects a significant number of miRNA-disease association entries from literature. Identified miRNAs were then searched for their differentially expressed mRNA targets using Cytoscape 55 through its app-CyTragetLinker.