Temporal Transcriptomic and Proteomic Landscapes of Deteriorating Pancreatic Islets in Type 2 Diabetic Rats

Type 2 diabetes (T2D) is a major public health issue characterized by pancreatic islet β-cell failure in the presence of insulin resistance. Accumulating evidence suggests that progressive deterioration of pancreatic β-cell function and gradual loss of β-cell mass could be the core events during T2D development, regardless of therapy status (1–4). Genome-wide association and sequencing studies have identified multiple risk variants for T2D, the majority of which appear to play a primary role in β-cell function rather than to affect insulin resistance, further highlighting the importance of β-cells in the pathogenesis of T2D (5).

T2D is a complex disease, and β-cell failure is likely caused by altered expression of many genes and their products. Therefore, the use of system-oriented strategies is critical to investigate the complex changes that occur in β-cells or pancreatic islets, which primarily comprise β-cells. Hence, large-scale and unbiased omics technologies, particularly microarray-based transcriptomics and mass spectrometry (MS)–based proteomics, have been used to analyze islets isolated from various T2D animal models and human cadaver donors to elucidate the mechanisms underlying β-cell failure (summarized in Supplementary Table 1). β-Cell failure during diabetes progression is a gradual process that undergoes various stages (6,7) wherein different molecule events occur in chronological order. However, current studies have focused primarily on single time points at relatively late stages of the disease, so mapping the order in which these events occur and distinguishing causal molecular events (leading to diabetes) from those that occur as a consequence of glucolipotoxicity associated with diabetic conditions are impossible. For this reason, prospective studies investigating the evolution of molecular events in islet β-cells at various stages of T2D are of interest.

The study of β-cells in humans with T2D often has been hindered by the limited accessibility of human islets and by ethical considerations. In this context, appropriate rodent models are essential for the identification of diabetic mechanisms (8). The Goto-Kakizaki (GK) rat, one of the best-characterized animal models of spontaneous T2D (9), shares many characteristics with human patients with diabetes (10). Similar to human T2D, the core cause underlying hyperglycemia in GK rats is β-cell failure (11,12).

In the current study, to understand the process of deteriorating pancreatic islets at the molecular level, we used RNA sequencing (RNA-seq) and tandem mass tag (TMT)–based quantitative proteomics technology to generate integrated transcriptomic and proteomic profiles of pancreatic islets in GK rats after the establishment of hyperglycemia (from 4 to 24 weeks). Subsequent bioinformatics analysis in a time course fashion revealed the chronological order of T2D-related molecular events during the deterioration of pancreatic islets. This large quantitative data set represents a valuable resource that provides a comprehensive picture of the mechanisms responsible for islet dysfunction and will allow us to identify potential interventions to prevent β-cell failure and deterioration in human T2D.

Brief descriptions of key experimental procedures are provided below. More details are provided in the Supplementary Data.

Founders of GK/Jcl diabetic rats were purchased from RIKEN BioResource Center (Ibaraki, Japan). All GK/Jcl diabetic rats and Wistar (WST) rats were maintained under specific pathogen-free conditions and were used between 4 and 24 weeks of age in accordance with the animal experimental guidelines set forth by the Institutional Animal Care and Use Committee of the Institute of Biophysics, Chinese Academy of Sciences.

Pancreatic islets from male GK and age-matched control WST rats were isolated through collagenase digestion. After separation on a Ficoll density gradient, the islets were handpicked in Hanks’ buffer under a dissection microscope.

Sixteen male GK rats were used in the following experiments. Littermates of GK rats were randomly divided into N-acetyl-L-cysteine (NAC) and control groups. Four-week-old rats were orally administered NAC 200 mg/kg of body weight (616-91-1; Sigma) or drinking water by gavage once a day for 12 weeks. Random blood glucose assay and glucose tolerance, insulin tolerance, and glucose-stimulated insulin secretion (GSIS) tests were performed as described in the Supplementary Data.

RNA-seq was performed by using a multiplex analysis of polyA-linked sequences (MAPS) approach as previously described (13).

Proteins extracted from isolated islets were digested and labeled by 6-plex TMT reagents (14) according to the instructions from the manufacture (Thermo Fisher Scientific). TMT-labeled peptide mixtures were equally pooled, separated by offline high pH reversed-phase chromatography, and repeatedly analyzed using a nano–liquid chromatography-tandem MS technique (Supplementary Fig. 1A). The raw MS data were processed with Proteome Discoverer 1.4 software. The peptide confidence value was set as 0.01. At the protein level, a precursor intensity fraction of 50% was selected as an optimal trade-off value for both identification and quantification (Supplementary Fig. 1B). A pseudocount representing relative protein abundance was calculated by using the TMT ratio and the normalized spectral abundance factor (15).

Both mRNA raw counts and protein pseudocounts were normalized by using the remove unwanted variation approach (Supplementary Fig. 1C) (16). Differentially expressed (DE) genes were assessed by ANOVA with a false discovery rate of <0.01. The Database for Annotation, Visualization and Integrated Discovery (DAVID) Web service API Perl Client (17,18) was used to perform gene ontology (GO) functional enrichment analysis (with a false discovery rate of <0.05). The k-means clustering algorithm was used to classify dynamic gene expression patterns. Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway enrichment analysis was carried out by using the Generally Applicable Gene Set Enrichment (GAGE) package in R software (with an adjusted P < 0.05) (19).

For this work, 92.4 gigabyte sequencing data were generated. All RNA-seq data were deposited in the National Center for Biotechnology Information Gene Expression Omnibus under accession number GSE 81811. The MS data were deposited in the ProteomeXchange Consortium through the PRIDE (Proteomics Identifications) (20) partner repository with the data set identifier PXD004709.

To investigate the global molecular dynamics of T2D islets, we analyzed the transcriptomes and proteomes of pancreatic islets isolated from male GK rats and age- and sex-matched control WST rats at five consecutive time points (weeks 4, 6, 8, 16, and 24) (Fig. 1A). Transcriptomes and proteomes of islets were measured by using the MAPS-based RNA-seq technique and TMT labeling–based proteomic method, respectively. Combined analysis of all samples yielded the identification of 15,101 mRNAs and 8,362 proteins, of which 7,395 overlapped (Fig. 1B). Furthermore, 13,866 mRNAs (minimal counts of 10 detected in at least three samples) and 5,631 proteins (identified in at least two biological replicates by peptides of precursor intensity fraction ≤50%) were considered a quantifiable data set, of which 5,015 overlapped (information for all identified genes is provided in Supplementary Table 2). We estimated relative protein abundance, which was noted as the protein pseudocount in this study, on the basis of TMT ratio and normalized spectral abundance factor. Compared with TMT ratios, protein pseudocounts resulted in the generation of more-reasonable clusters that matched the experiments (Supplementary Fig. 2A). Protein pseudocounts positively correlated with mRNA counts, with a mean Spearman correlation coefficient of 0.37 (Supplementary Fig. 2B), which is similar to reported values in other biological systems (21,22). These results demonstrate that our method of estimating protein pseudocounts was reasonable and unbiased.

We first carried out an unsupervised hierarchical clustering analysis of the transcriptomes and proteomes from GK and WST islets (Fig. 1C). At the mRNA level, WST islets clustered together and clearly separated from GK islets as expected, representing different pathophysiological states of islets in control rats versus diseased rats. Overall, the transcriptomes of GK and WST rats demonstrated appropriate clustering at different time points, representing the developmental stages of islets. Of note, GK islets at 4, 6, and 8 weeks formed one branch that was distinct from the branch formed at 16 and 24 weeks, likely reflecting two different diabetic stages of islets in GK rats. Upon proteome clustering examination, GK islets at 4 and 6 weeks were separated from other GK islets and instead clustered with WST islets, suggesting minimal changes in protein expression at the early stages of T2D, even when there are significant changes at the mRNA level.

We next performed a principal component (PC) analysis to investigate transcriptome and proteome dynamics over time (Fig. 1D). PC2 primarily reflected age-related developmental changes, whereas PC1 represented the differences between normal WST and diabetic GK rats. Similar to the above results obtained from unsupervised hierarchical clustering, PC1 highlighted two notable diabetic stages in GK islets. At the mRNA level, islets at 4, 6, and 8 weeks clustered as one stage, and islets at 16 and 24 weeks represented another stage. However, at the protein level, islets at 4 and 6 weeks clustered as one stage, and the remaining time points clustered as a separate stage. GK islets progressively develop into disorganized structures exhibiting pronounced fibrosis separating strands of endocrine cells (23). Of note, these changes were not present or rare in islets at 4–6 weeks but became prominent at older ages (8–24 weeks) (Supplementary Fig. 2C), correlating well with our proteome-defined two stages (Fig. 1C and D). Taken together, global profiling at the mRNA and protein level roughly characterizes the deterioration of islets in GK rats from 4 to 24 weeks of age into two stages: an early stage at 4–6 weeks and a late stage at 8–24 weeks, with a turning point at ∼8 weeks.

DE genes between GK and WST rats at each time point were first assessed by one-way ANOVA (P < 0.05) (Supplementary Table 2). The results revealed a remarkable increase of DE mRNAs over time, whereas only a mild increase was found in the number of DE proteins, suggesting a more dynamic regulation of gene expression at the mRNA level than at the protein level during the development of GK diabetes (Supplementary Fig. 3A). When we compared the fold changes in mRNA and protein levels, the Pearson correlation value (0.11) was relatively poor at 4 weeks but increased to 0.23–0.28 at later time points (Supplementary Fig. 3B).

From our time-resolved expression data set, we analyzed the temporal significance of gene expression changes by using two-way ANOVA, which considered weeks (five different time points) and rats (GK vs. WST) as the two statistical factors. In total, we identified 5,551 DE genes, including 3,910 mRNAs and 2,387 proteins (Supplementary Table 2), of which 746 were identified at both the mRNA and the protein level. The correlations among these 746 DE genes varied from anticorrelation to full accordance, with an average Pearson coefficient of 0.39 (Fig. 2A and Supplementary Table 3), which is larger than the correlation coefficients at individual time points. DAVID functional clustering analysis indicated carbon metabolism and ribosome enrichment among genes with concordant mRNA and protein levels (Fig. 2B), potentially representing the set of genes exhibiting stable mRNA and protein expression (24). In contrast, no significant functional clustering was identified for negatively correlated DE genes.

For the confirmation of expression data, we randomly selected six DE genes to measure their mRNA expression by quantitative RT-PCR (qRT-PCR) in independent GK/WST islet samples. The expression patterns of these genes measured by RNA-seq and qRT-PCR were similar (Supplementary Fig. 3C). Moreover, by comparing our GK/WST data set with the published data set of islets from individuals with and without T2D (Supplementary Table 1) (10), we found that on average, 68.9% of DE genes were consistent between GK rats and humans with diabetes (Supplementary Fig. 3D), indicating high relevance of the current study in GK rat to human islets in T2D.

The primary purpose for the generation of the time course data set in this study is to reveal the temporal properties of biological pathways relevant to the development of GK diabetes at the system level. Therefore, we performed time course pattern analysis for all DE genes by using the k-means clustering method and successfully identified 12 mRNA expression patterns (m1–m12) and 9 protein expression patterns (p1–p9) (Fig. 2C). To explore the biological functions of these expression patterns, we carried out a DAVID analysis and organized the identified networks with enriched functions by using EnrichmentMap software (Supplementary Table 4) (25). On the basis of the time-course expression patterns of clustered genes, we classified the biological events into the following categories:

1. Constant up (m2, p8) and down (p2, m3, m10) genes, which were either up- or downregulated at all the time points. The constant up genes were highly enriched for such functions as cell redox homeostasis, translation elongation, cytoskeleton organization, antiapoptosis, and so forth. In contrast, the constant down genes were mainly associated with mitochondrion, metabolism, lysosome, protein transport, and so forth.

2. Up early genes (p1, p6, m11) that were upregulated at 4–8 weeks, which include those participating in the glucose metabolism and innate immune response. In addition, many proteasome proteins were dramatically upregulated at 4 weeks.

3. Up late genes (p9, m8, p7, m1, m12) that were upregulated at the late stage of 8–24 weeks. These genes are highly enriched for cell adhesion and cytoskeleton organization at the protein level, probably associating with the development of islets fibrosis, whereas at the mRNA level, apoptosis was significantly upregulated at 24 weeks.

4. Down early genes (p3, m6) that were downregulated at 4–8 weeks. The most noticeable feature is the downregulation of oxidative phosphorylation (OXPHOS) and tricarboxylic acid (TCA) cycle at the protein level, probably suggesting the insufficient energy supply and oxidative stress in GK islets at the early stage (26–28). At the mRNA level, the genes associated with cell cycle and nuclear lumen were highly enriched, likely contributing to the loss of β-cells in GK islets.

5. Down late genes (p5, m9, p4, p7, m5, m4) that were downregulated at 8–16 weeks. The GO functions of insulin secretion, lysosome, and secretory granule were representatively enriched.

Such a time course clustering analysis from genes to biological functions suggests that the progression of GK islet deterioration was stage based and dynamically regulated. Furthermore, because genes exhibiting similar expression patterns generally share functional relationships, clustering analyses also allow the prediction of genes that share similar temporal expression patterns, with previously validated diabetes-related genes as potential new candidates for further investigation.

To gain a deeper understanding of temporal pathway sequences during the deterioration of GK islets, we performed GAGE analysis with the quantitative transcriptomic and proteomic data sets and identified 161 KEGG pathways significantly enriched for at least one time point (Benjamini-Hochberg–adjusted P < 0.05) (Supplementary Table 5). Consistent with the above gene expression clustering analysis results, these KEGG pathways also roughly comprised down- and upregulation-dominated temporal classes (Fig. 3). For example, insulin secretion, and SNARE interactions in vesicular transport were downregulated at the late stage with the same temporal pattern (representative genes illustrated in Supplementary Fig. 4). These two pathways were partially downregulated at 4 and 6 weeks with no significant change in statistics and became significantly downregulated after 8 weeks, clearly demonstrating the dynamic features of various pathways involved in insulin secretion during the progression of T2D. We also noticed that the pathway of glycolysis/gluconeogenesis was gradually upregulated since 6 weeks, probably indicating that in GK islets, the anaerobic metabolism is the dominant approach for supplying the energy because of the defect in OXPHOS and TCA cycle.

In addition, we found that signaling pathways associated with inflammation were upregulated at the mRNA level, including the nucleotide oligomerization domain (NOD)–like receptor, tumor necrosis factor (TNF), and nuclear factor-κB (NF-κB) signaling pathways. This finding is consistent with previous reports that islet inflammation plays an important role in the pathogenesis of GK diabetes (26,29) and provides more comprehensive details of pathway dynamics at both the mRNA and the protein level.

We identified 311 DE genes as mitochondria related by using GO terms for cellular components that contain the keywords mitochondrion or mitochondrial and divided them into four groups on the basis of the unsupervised hierarchical clustering analysis of their temporal profiles (Supplementary Table 6). Further manual annotation revealed the details regarding mitochondrial dysfunction during the progression of T2D (Fig. 4). We found that OXPHOS complexes, mitochondrial ribosome proteins, translocase outer/inner membrane complex (TOM/TIM), and some metabolite transporters were downregulated early at the protein level, and conversely, most of corresponding mRNAs were upregulated at 4 and 6 weeks, indicating transcription compensation at the early stage of T2D. However, this compensation ability was eventually lost at the later stage during β-cell deterioration. In addition, several proteins responsible for protein assembly and quality control, mitochondrial biogenesis, and mitochondrial DNA transcription were downregulated, such as Hspd1, Grpel1, Lonp1, Letm1 and Pmpcb, Tfam, Mtfr1l, Mfn2, the transfer RNA ligases (Rars2, Nars2, Dars2, and Vars2), and 12 mitochondrial RNAs (Supplementary Fig. 5A). Taken together, mitochondrial dysfunction was considered one of earliest pathogenic events in islets of GK rats.

To gain a deeper understanding about how the metabolism in GK islets changed with the development of diabetes, we mapped our quantitative omics data to KEGG metabolic pathways (Fig. 5). The most notable change of metabolism in GK islets was the upregulation of glycolysis metabolism (26,27) and the downregulation of the TCA cycle, OXPHOS, and fatty acid metabolism, which suggests that the primary metabolism of GK islets switches from aerobic metabolism to anaerobic metabolism (the so-called Warburg-like effect). Furthermore, we found that besides FAD-dependent glycerol-3-phosphate dehydrogenase (Gpd2) reported previously (30), Got1, Got2, Mdh1, and Mdh2 were also downregulated, representing a comprehensive picture of defective malate-aspartate shuttle and glycerol phosphate shuttle in GK islets (Supplementary Fig. 5B). Such a defect could cause an increase of the cytoplasmic NADH/NAD⁺ ratio, enhance the formation of lactate from pyruvate, and further disrupt the link between cytosolic glycolysis and mitochondrial metabolism.

In addition, some genes participating in amino acid metabolism were also downregulated, partially contributing to a defect in GSIS (31–33). However, several enzymes involved in glutathione metabolism were upregulated at both the mRNA and the protein level in GK islets, including Gclc, Ggct, Ggt1, Gsta3, and Gsto1, which may reflect an adaptive mechanism to combat increased reactive oxygen species (ROS) in GK islets (34,35).

Insufficient insulin secretion is caused by either impaired GSIS or reduced β-cell mass. β-Cell mass is regulated by a balance among β-cell neogenesis, proliferation, and apoptosis. Consistent with previous studies (26,36), the current data show that several genes functionally associated with neogenesis are downregulated in GK islets, including Pdx1, Nkx2-2, Nkx6-1, Mafa, and Fev, at the mRNA and/or protein level (Fig. 6). Moreover, several genes specific for α-cells and pancreatic polypeptide cells (i.e., Arx, Isl1, Pax6, Pou3f4, Ppy) were downregulated at the early stage. Certain genes required for the differentiation of pancreatic progenitors into endocrine progenitors (i.e., Foxa1, Gata6, Sox9, Onecut1) were significantly upregulated at the late stage of GK diabetes.

Low levels of β-cell proliferation in GK rats constitute another factor contributing to decreased β-cell mass (11,26,36). In the current data set, 38 genes among cluster m6, including Aurkb, Ccna2, and Kifc1, were associated with cell cycle and nuclear lumen and exhibited similar expression patterns at the early stage (Supplementary Fig. 6A). These genes gradually decreased over time in WST rats, reflecting an aging-dependent reduction in proliferation over time (38). In contrast, all these genes were dramatically downregulated at 4 weeks in GK islets and then progressively decreased to even lower levels at 24 weeks (Supplementary Fig. 6B). Consistently, Ki67-positive β-cells were significantly decreased in GK islets at 4, 6, and 8 weeks (Supplementary Fig. 6C), suggesting reduced β-cell proliferation (11,12,26). Conversely, the apoptosis pathway was only elevated at 24 weeks at the mRNA level, implying that apoptosis was not responsible for β-cell loss during the early phase of the disease.

Chronic inflammation in GK islets has been demonstrated and considered as a pathophysiological contributor in T2D (29,39). In the current study, temporal expression of proinflammatory cytokines revealed two distinct stages: early priming and late amplification. During the priming stage, interleukin-1β (IL-1β) and IL-6 were only elevated to low levels between 4 and 8 weeks followed by a rapid increase to much higher levels (84-fold increase in GK at 24 weeks for IL-6) during the late amplification stage (between 16 and 24 weeks) (Fig. 7). Because immune cell infiltration was hardly detectable in GK islets before 8 weeks (39), the early priming stage is likely induced by metabolic dysfunction in GK islets. The identity of specific sensors that are triggered to produce this priming of inflammation is not fully understood.

The NLRP3 (NLR family, pyrin domain containing 3) inflammasome activates both NF-κB (to induce pro-IL-1β production) and caspase-1 (to process pro-IL-1β into its mature active form). ASC (apoptosis-associated speck-like protein containing a CARD, also called PYCARD), which interacts with NLRP3 during inflammasome assembly, was upregulated in GK islets, indicating possible activation of the NLRP3 inflammasome during the initiation of sterile inflammation. The NLRP3 inflammasome is also activated by TXNIP (thioredoxin-interacting protein), a key node linking glucotoxicity and endoplasmic reticulum stress to NLRP3 inflammasome activation (40). We observed that TXNIP was gradually upregulated in GK islets after 6 weeks.

Oxidative stress is a pathogenic factor caused by chronic hyperglycemia in GK pancreatic islets (41). However, sources of ROS in T2D have not been clearly defined (42). The current data suggest multiple sources of ROS accumulation in GK islets at the early stage, including dysfunctional mitochondria, nitric oxide synthase (Nos2/iNOS), NADPH oxidase (Nox4), cyclo-oxygenase (Ptgs1/Cox1 and Ptgs2/Cox2), and cytochrome P450 mono-oxygenase (Cyp2s1, Cyp7b1, and Cyp4f5) (Fig. 7). We also observed increased expression of several antioxidants, such as Sod1, Prdx4, and Gpx2 (Fig. 7), in GK islets, consistent with previous reports (27,35,43).

To validate the involvement of ROS in the pathogenesis of T2D, we treated 4-week-old GK rats with an antioxidant, NAC, for 10 weeks. NAC treatment significantly reduced random blood glucose of GK-NAC rats (n = 8, NAC 200 mg/kg) compared with that of age-matched GK-control rats (n = 8, sham) (Supplementary Fig. 7A–D), and GK-NAC rats exhibited greater tolerance to high glucose in the oral glucose tolerance test. No significant differences in insulin sensitivity between GK-NAC and GK-control rats were observed. The islets in GK-NAC groups had higher GSIS than those in GK-control groups (Supplementary Fig. 7E). Furthermore, we measured the mRNA expressions in islets between GK-NAC and GK-control by using RNA-seq. The results show that the insulin secretion pathway was upregulated, whereas ROS and inflammation-related genes (Tnf, Ggt1, Pycard, Cxcl1, etc.) and signaling pathways (NOD-like receptor signaling pathway, TNF signaling pathway, metabolism of xenobiotics by cytochrome P450, etc.) were downregulated in GK islets after NAC treatment (Supplementary Fig. 7F). Thus, the neutralization of increased ROS by NAC ameliorates its impact on GSIS and inflammation and protects GK β-cell function.

We carried out a large-scale analysis of gene and protein dynamics in pancreatic islets of GK rats at various stages of T2D. Combined transcriptome and proteome analysis revealed sufficient depth of coverage and quantitative accuracy to generate functional portraits of healthy and diseased pancreatic islets with unprecedented detail. Many DE genes and proteins identified previously were confirmed in this study, such as those associated with OXPHOS, mitochondrial function, metabolism, insulin secretion, oxidative stress, and inflammation, thus validating the analytic methods used here. More importantly, the construction of time course–based gene expression and protein profiles allowed us to identify the chronological order of biological events contributing to the pathogenesis of T2D.

The data suggest two early events that likely contribute to T2D in GK islets: a reduction in β-cell mass and a shift in metabolism. Many transcription factors required for the specification of endocrine cells were downregulated early in GK islets, whereas those required for trunk and exocrine cells were not altered at 4 weeks but increased at later time points (Fig. 6). Indeed, GK rats from the Paris colony exhibit a significant reduction in β-cell mass at the fetal stage that precedes the onset of hyperglycemia at ∼4 weeks after birth (11), similar to our GK colony. Besides neogenesis, defective proliferation also causes a reduction in β-cell mass, as has been suggested for GK rats from the Paris colony (11,26). In our GK rats, many genes required for the cell cycle and proliferation were significantly downregulated at the early stage (4–6 weeks) (Supplementary Fig. 6C and D), suggesting an early defect in proliferation.

Another notable feature of the current data is the observation of an early shift in metabolism. As early as 4 weeks, the primary metabolism in the islets of GK rats switched from aerobic metabolism to anaerobic metabolism (Warburg-like effect). Although metabolism switch has been proposed in previous research (34), our time course–based quantitative data allowed us to detect this metabolic shift as an early event of T2D in GK rats and permitted us to gain a deeper insight into the mechanism underlying this shift. Insufficient reduction in β-cell mass alone may not necessarily cause T2D because autopsy studies of patients with T2D have revealed an ∼50% decrease in β-cell mass compared with BMI-matched control patients (44,45). Furthermore, a reduction in β-cell mass of ∼50% is required for dogs and rats to develop diabetes (46). This metabolic shift is likely caused by mitochondrial dysfunction. Many proteins associated with the OXPHOS system, metabolite transporters, and the TCA cycle were downregulated starting at 4 weeks (Fig. 3). Thus, the mitochondrial limitation of glucose oxidation in GK islets occurred during the early stage. To compensate for this energy deficiency, GK islets improved the rate of glycolysis and upregulated the expression of Ldha, which converted pyruvate to lactate and further disrupted the link between cytosolic glycolysis and mitochondrial metabolism (Fig. 4). Anaerobic glucose metabolism with NADH accumulation in the β-cell of mitochondrial diabetes caused by ethidium bromide treatment can impair the transcription of mitochondria DNA, halt the TCA cycle, and affect GSIS (47). Therefore, this Warburg-like metabolic shift may contribute to the early impairment of GSIS in GK β-cells because GSIS requires the production of both a triggering signal (ATP) and amplifying signals (i.e., cAMP, short-chain acyl-CoA compounds, NADPH) produced during aerobic metabolism (48).

Despite the early reduction in β-cell mass, insulin and the enzymes required for proinsulin processing were comparable or even higher at both the mRNA and the protein level at 4–6 weeks, suggesting compensation during the early phase in response to hyperglycemia. Absolute insulin insufficiency only occurred during the late stage of T2D. Consistently, we also observed compensation for OXPHOS complexes and mitochondrial machinery at the mRNA level at 4–8 weeks, despite reduced protein levels. Thus, the progression of T2D in GK rats from 4 to 24 weeks can be divided into two stages: compensation (4–6 weeks) and compensation failure (8–24 weeks). Further data mining to examine temporal expression patterns will help to elucidate the mechanisms underlying compensation and decompensation.

Mitochondrial dysfunction and the Warburg effect generate greater ROS invasion, which in turn induces chronic low-grade inflammation (49). Of note, proinflammatory cytokines also exhibit two distinct stages (Fig. 7): a priming stage at 4–6 weeks and an amplification stage after 8 weeks. Given that no inflammatory cell infiltration was observed in GK islets before 8 weeks (39), the priming stage of inflammation was likely induced by intracellular signals generated by metabolic stress, such as the presence of increased ROS or free fatty acids. However, the second stage of inflammation amplification may be induced by a complex combination of intracellular and extracellular (i.e., macrophage infiltration) inducers. The current data provide clues to unravel the mechanism underlying the initiation and amplification of sterile inflammation; understanding this mechanism is necessary to developing novel anti-inflammation therapies to treat T2D. Islet inflammation is undoubtedly an early event during T2D pathogenesis, but it is not likely a causal event because IL-1β and TXNIP were not significantly expressed at 4 weeks, although they gradually increased during later stages of the disease.

In summary, the data reveal two stages during the progression of T2D in GK islets. The early stage (4–6 weeks) is characterized by anaerobic glycolysis, inflammation priming, and compensation for insulin synthesis, whereas the late stage (8–24 weeks) is characterized by inflammation amplification and compensation failure (Supplementary Fig. 8). We did not observe significant apoptosis during the early stage. The apoptosis pathway was only significantly elevated at 24 weeks at the mRNA level. The time course transcriptome and proteome data sets for GK rat islets depict a comprehensive landscape of dynamic changes in gene expression at various stages of diabetes, representing a valuable resource for the research community to further explore the molecular etiology and progression of diabetes. In-depth exploration of this resource will aid in the discovery of potential diagnostic and therapeutic targets for human T2D.

Acknowledgments. The authors thank the staff of the Institute of Biophysics Core Facilities, in particular, Yan Teng for technical support with confocal imaging, Dr. Jifeng Wang for MS operation, and Zhen Fan and Xiaowei Chen for RNA-seq design and data collection.

Funding. This work was supported by grants from the National Key Basic Research Project of China (2014CB910503 and 2015CB910303), the National Key Research and Development Program of China (2016YFC0903301), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA12030101), and the National Science Foundation of China (31421002, 31400703, and 31400658).

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. J.H. performed the proteomic experiments and MS data analysis. J.H. and W.Z. prepared the figures. J.H., W.Z., and D.Z. carried out the bioinformatics analyses. J.H., Y.W., and T.X. wrote the manuscript with help from the other authors. Z.L. performed the RNA-seq experiments and immunohistochemistry imaging. Q.H. carried out rat breeding. L.L. and L.W. performed the qRT-PCR and GSIS experiments. Y.W. carried out the islet preparation and animal experiments. Y.W. and T.X. conceived the project. All authors read the manuscript and discussed the interpretation of results. T.X. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.