Temporal Transcriptome Analysis Reveals Dynamic Gene Expression Patterns Driving β-Cell Maturation

Newly differentiated pancreatic β cells lack proper insulin secretion profiles of mature functional β cells. The global gene expression differences between paired immature and mature β cells have been studied, but the dynamics of transcriptional events, correlating with temporal development of glucose-stimulated insulin secretion (GSIS), remain to be fully defined. This aspect is important to identify which genes and pathways are necessary for β-cell development or for maturation, as defective insulin secretion is linked with diseases such as diabetes. In this study, we assayed through RNA sequencing the global gene expression across six β-cell developmental stages in mice, spanning from β-cell progenitor to mature β cells. A computational pipeline then selected genes differentially expressed with respect to progenitors and clustered them into groups with distinct temporal patterns associated with biological functions and pathways. These patterns were finally correlated with experimental GSIS, calcium influx, and insulin granule formation data. Gene expression temporal profiling revealed the timing of important biological processes across β-cell maturation, such as the deregulation of β-cell developmental pathways and the activation of molecular machineries for vesicle biosynthesis and transport, signal transduction of transmembrane receptors, and glucose-induced Ca²⁺ influx, which were established over a week before β-cell maturation completes. In particular, β cells developed robust insulin secretion at high glucose several days after birth, coincident with the establishment of glucose-induced calcium influx. Yet the neonatal β cells displayed high basal insulin secretion, which decreased to the low levels found in mature β cells only a week later. Different genes associated with calcium-mediated processes, whose alterations are linked with insulin resistance and deregulation of glucose homeostasis, showed increased expression across β-cell stages, in accordance with the temporal acquisition of proper GSIS. Our temporal gene expression pattern analysis provided a comprehensive database of the underlying molecular components and biological mechanisms driving β-cell maturation at different temporal stages, which are fundamental for better control of the in vitro production of functional β cells from human embryonic stem/induced pluripotent cell for transplantation-based type 1 diabetes therapy.

Introduction

Pancreatic β cells are functionally defined by their capacity for insulin secretion, stimulated by glucose and other nutrients. Loss of functional pancreatic β cells is the primary cause of diabetes, and researchers have intensively studied β-cell development for the last two decades to generate new therapeutic approaches. Type 1 diabetes results from autoimmune destruction of β cells in the pancreatic islet, whereas the more common type 2 diabetes results from peripheral tissue insulin resistance and β-cell dysfunction. Diabetic patients, particularly those suffering from type 1 diabetes, could potentially be cured through transplantation of new β cells. To this end, several protocols have allowed the production of glucose responsive β-like cells from human embryonic stem/induced pluripotent cells (Kushner et al., 2014; Schiesser and Wells, 2014; Tabar and Studer, 2014). These β-like cells show gene expression, ultrastructural characteristics, and glucose responsiveness, both in vitro and in vivo, which closely resembling the features of β cells found in pancreatic islets (Pagliuca et al., 2014). However, the production of these cells is limited as the final cell population has about 30–60% β-like cells, and many of the remaining cells are relatively uncharacterized and can be undifferentiated progenitors or other types of unwanted cells (Shahjalal et al., 2018). This low efficiency is partly due to our lack of understanding the signaling pathways that direct β-cell maturation (Kieffer, 2016).

Newly made insulin-expressing cells, in both the human and rodent fetus, are not mature (pre-β or immature β). They secrete two to five times more insulin than adult β cells with basal glucose (<5.6 mM) while lacking robust glucose-stimulated insulin secretion (GSIS) under stimulating (>10 mM) glucose (Adam et al., 1969; Rorsman et al., 1989; Hellerstrom and Swenne, 1991; Bliss and Sharp, 1994; Rozzo et al., 2009). A maturation process converts pre-β cells into mature β cells with low basal insulin secretion but high GSIS. Several molecular mechanisms can promote β-cell maturation: insulin biosynthesis and vesicle packaging are necessary for insulin secretion (Gu et al., 2010; Blum et al., 2012; Goodyer et al., 2012); glucose influx into β cells, glycolysis, and oxidative phosphorylation lead to ATP production, which represses ATP-sensitive potassium channels to induce Ca²⁺ influx and trigger insulin secretion (Hou et al., 2009; Henquin, 2011); intercellular communication controls the coordinated and pulsatile nature of insulin secretion via GAP junctions (Benninger and Piston, 2014) or heterotypic protein interactions (Konstantinova et al., 2007); several nutritional and neural signals are established to control the dose of secretion to properly regulate blood sugar for physiological demands, thus requiring the production and subcellular localization of neural transmitter receptors and effector channels (Osundiji and Evans, 2013). Different transcriptional factors and signaling molecules, including MafA, NeuroD, and calcineurin (Nishimura et al., 2006; Gu et al., 2010; Goodyer et al., 2012), control and promote the maturation processes. Similarly, gene regulatory mechanisms for proper expression of metabolic genes, such as DNA methylases, regulate maturation by controlling proper glucose metabolism (Dhawan et al., 2015). Despite these published studies, many questions on maturation remain unresolved. Specifically, the reported stage of maturation varies from 2 days to 2 weeks after birth in rodents (Nishimura et al., 2006; Rozzo et al., 2009; Blum et al., 2012), and it is unknown which mechanism(s) represent(s) the limiting step for maturation. To this end, several studies compared the gene expression profiles between mature and immature β cells (Nielsen et al., 2004; Jermendy et al., 2011; Blum et al., 2012; Benitez et al., 2014; Hrvatin et al., 2014; Dhawan et al., 2015). Yet, these comparisons did not consider the dynamics from progenitor to mature β cells, which are necessary to distinguish genes associated with β-cell differentiation (production of insulin⁺ cells) and/or maturation (gaining glucose response). Only few studies monitored gene expression of β-cell maturation at multiple stages so far and recently also at single-cell level, using computational methods able to provide a pseudotemporal ordering of the cells (Qiu et al., 2017). However, in these studies, the altered gene expression levels across the stages are analyzed by computational methods for differential expression such as DEseq2 (Love et al., 2014), comparing each time point independently and without considering the temporal profile of each gene. In addition, the biological interpretation of the obtained lists of differentially expressed genes is usually performed through enrichment techniques, which are applied independently from the gene selection step, thus generating potential false positives/negatives. Finally, the expression profiling of the selected genes is typically displayed through heatmaps, mostly dichotomized into immature and mature cells, without characterizing clusters of temporal expression profiles representing the dynamics of β-cell development and maturation through all the stages.

Here, we examined the dynamics of gene expression across six key stages of β-cell maturation, starting from precommitted endocrine progenitors to adult functional β cells. The aim of the study was to determine temporal patterns (TPs) of functional groups of genes that promote newly born, nonfunctional β cells to become functional glucose-responsive β cells. With respect to the previous analyses, a computational pipeline recently published, named FunPat, was applied to build a comprehensive map of the temporal evolution of functional processes and genome-wide candidate markers. Specifically, FunPat combines gene selection, clustering of temporal expression profiles, and functional annotation into a single framework, and it has shown high precision and recall in detecting the correct temporal expression patterns (Sanavia et al., 2015). The resulting dynamic gene expression profiles were then correlated with the temporal development of β-cell properties associated with GSIS, including insulin vesicle biosynthesis, glucose-induced calcium influx, and insulin secretion. Finally, we also provided a comprehensive database describing the biological processes and pathways characterizing β-cell maturation across time. Among them, several established gene products involved in calcium signaling, whose deregulation critically affects homeostasis in insulin-secreting β cells, showed significant temporal correlation between their enhanced expression and maturation, both highlighting potential targets for diseases such as diabetes and providing useful insights for in vitro derivation of β cells to develop new therapies.

Materials and Methods

Animal Breeding and Use

Mouse usage followed the procedures specified in protocols M/10/168 and M/11/181 approved by the Vanderbilt Institutional Animal Care and Use Committee. Mip^eGFP, Ngn3^eGFP, and Rip^mCherry mice were from Hara (Hara et al., 2003), Kaestner (Lee et al., 2002), or reported (Zhu et al., 2015), respectively. The stock mice were kept in mostly Bl6/CBA and 129Sv/Ev mixed background with intercrosses (P0). When heterozygous mice were needed, the stock breeders were directly crossed with CD1 mice to produce P1 progeny for pancreatic tissue collection (Charles River). For collecting Ngn3^eGFP/eGFP homozygous embryos, intercrosses between P1 mice were utilized. This crossing scheme allows collection of mice with similar genetic background.

Antibody Staining

Antibody staining followed routine protocols. Rabbit antiamylase (Sigma–Aldrich) and guinea pig anti-insulin (Dako) were used at 1:1,000 dilution. Cy5-conjugated donkey anti–guinea pig and fluorescein isothiocyanate (FITC)–conjugated donkey anti-rabbit (Jackson Immunology) were used at 1:2,000. Hand-picked islets were used for whole-mount staining and imaged with confocal microscopy.

Cell Sorting, RNA Extraction, Real-Time Polymerase Chain Reaction, and RNA Sequencing

To examine the gene expression dynamics along the maturation steps, we used four β-cell populations: P1 (day 1 after birth) and P4 β cells are immature cells (Blum et al., 2012); P12 and P60 β cells represent newly matured and fully functional β cells, respectively (Nishimura et al., 2006; Blum et al., 2012; Dhawan et al., 2015). Two progenitor pools transcribing Ngn3 at embryonic stage E15.5 without (Ngn3^eGFP/eGFP) or with Ngn3 protein (Ngn3^eGFP/+) were included to establish a baseline of gene expression for differentiation. A portion of the progenitors will become β cells (Gu et al., 2002).

EGFP⁺ and mCherry⁺ cells were sorted using Aria III (BD). Ngn3^eGFP pancreata were dissociated with trypsin (Gu et al., 2004) and used for the embryonic stages. For postnatal stages, Mip^EGFP and Rip^mCherry islets were first hand-picked, allowed recovering for 2 h in RPMI media [Life Technologies, 10% fetal bovine serum (FBS)], quickly washed with Ca-Mg free Hanks balanced salt solution (HBSS) (Cellgro) once, and dissociated with trypsin. For hand-picking, P1, P4, and P12 pancreata were directly digested with freshly made 1 mg/ml collagenase IV (Sigma–Aldrich) in HBSS. Perfusion was used for P60 pancreata. After digestion (<10 min), pancreatic clusters were quickly washed with RPMI. All solutions/media used throughout islet isolation have non-stimulating glucose at 5.6 mM. RNA was prepared with TRIzol (Life Technologies) and a DNA-free RNA^TM kit (Zymo Research). Two hundred monograms total RNAs with RINs greater than 8 (assessed using Agilent Bioanalyzer 2100) were sequenced, with three biological replicates, following Illumina protocols on HiSeq-2000. Real-time polymerase chain reaction (RT-PCR) utilized SYBR-Green (Promega), following manufacturer’s procedures.

RNA Sequencing Data Preprocessing and Function-Based Pattern Analysis

Raw reads were aligned to the mouse genome (mm10) and transcriptome using STAR version 2.3.0e (Dobin et al., 2013). Data were normalized and quantified with PORT pipeline to determine the relative expression level of each gene¹. FunPat was used to select differentially expressed genes that share functional annotation and common dynamic expression profiles (Sanavia et al., 2015). Genes were first filtered using the bounded-area method (Di Camillo et al., 2007), which calculates, for each gene, the area of the region bounded by the time series expression profile and a baseline, set at the expression level in E15.5 Ngn3^eGFP/eGFP cells as they have no active β cell–specific genes. A p-value was assigned to each gene by evaluating the significance of its bounded-area against a null hypothesis distribution described by a log-normal estimated comparing the biological replicates and characterized by mean and standard deviation equal to 1.09 and 0.44, respectively. Applying a Bonferroni correction to the resulting p-values, we named two sets of genes, seeds and candidates, selected considering adjusted and unadjusted p-values below 1%, respectively.

Both seeds and candidates underwent a function-based clustering approach (Di Camillo et al., 2012; Sanavia et al., 2015), which searches for TPs by performing a linear model-based clustering on the expression profiles, after subtracting the baseline, of groups of genes annotated to the same functional term, e.g., a Gene Ontology (GO) term (Ashburner et al., 2000) or a pathway (Kamburov et al., 2011). Each identified TP contains at least a seed, and it represents the mean differential expression across time with respect to the baseline, i.e., expression at E15.5 Ngn3^eGFP/eGFP. To lower the percentage of false negatives from the bounded-area method while preserving the false discovery rate, the algorithm selects a gene as differentially expressed if it is a seed or if it is a candidate that belongs to a TP and therefore shares the same biological annotation with at least a seed (Sanavia et al., 2015). Intuitively, all the genes associated with the same TP are likely to be differentially expressed as they are highly correlated to the same temporal profile and, as the clustering is functional term-specific, to a common biological function.

Functional Annotations

In order to achieve a comprehensive coverage of the functional terms and gene annotations currently available, both functional MGI (Mouse Genome Informatics) annotations to 15,939 GO terms from all categories (Biological Process, Molecular Function and Cellular Component)² and 1,628 pathways from Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome derived from ConsensusPathDB (Kamburov et al., 2011)³ were considered. While ConsensusPathDB already provides pathway annotations with reduced redundancy by mapping physical entities from different source databases to each other, in order to address the high information redundancy affecting GO annotation, FunPat exploits the hierarchical structure of GO database to search the TPs. Specifically, the method first starts from the most specific terms, represented by the leaf nodes in the ontology, and then it removes the differentially expressed genes associated with a significant TP from the annotations of the ancestor nodes, representing more general concepts (Sanavia et al., 2015). Finally, to summarize the results, the GO terms and pathways selected by FunPat were grouped into (1) functional categories related to common ancestor terms and (2) biological processes potentially related to GSIS according to MacDonald et al. (2005) and Benitez et al. (2012) and a manually-defined list of reference GO ancestor terms. Pathways were linked to the GO terms with the most similar meaning and associated with the corresponding common ancestor term. Pathways specifically related to diseases were grouped into a separate functional category named “Disease.”

Time-Dependent Organization of the Temporal Patterns

To better interpret the main dynamics characterizing the biological mechanisms involved in β-cell maturation, the identified TPs were sorted according to the most representative maturation stage. Specifically, based on the temporal progression x = {x₁ = E15.5 Ngn3^eGFP/eGFP, x₂ = E15.5 Ngn3^eGFP/+, x₃ = P1, x₄ = P4, x₅ = P12, x₆ = P60} and the TP = {TP(x₁),..,TP(x₆)}, one or more time break(s) was(were) assigned to each TP. A specific stage x_i was identified as a time break if, observing the pairs {x_i–1, x_i} and {x_i, x_i+1}:

$\frac{\left|\text{TP}(x_i)-\text{TP}(x_{i-1})\right|}{\overline{\text{TP}}}\cup \frac{\left|\text{TP}\left(x_{i+1}\right)-\text{TP}\left(x_i\right)\right|}{\overline{\text{TP}}}>0.2$

i.e., the TP showed in at least one of the pairs an increased or decreased level with a change over 20% with respect to the mean expression of the corresponding pattern. Time breaks were searched starting from E15.5 Ngn3^eGFP/+, whereas for the last stage (P60) the variation with respect to the time breaks identified in the previous stages was also evaluated. If more than one time break was assigned to the same TP, the break corresponding to the highest expression difference with respect to the baseline was selected, creating at the end five time-dependent groups of TPs. The TPs associated with the same time breaks were then summarized into main patterns (MPs) by applying the linear model-based clustering approach used in the FunPat pipeline. In this way, while TPs represent clusters of genes, the MPs are representative of clusters of functional terms, increasing the interpretability of the results.

Biological Interpretation of the Main Patterns

Each MP was classified into “positive” or “negative” according to the corresponding sign at its most representative time break, e.g., an MP = {MP(x₁),..,MP(x₆)} having the most representative time break at x₅ (i.e., P12) is defined positive if MP(x₅) > 0. Fisher’s exact test was applied to identify the statistically enriched biological functions, i.e., the functional categories or GSIS processes, in two cases: (1) genes belonging only to positive or negative MPs, regardless of timing and (2) focusing on positive or negative MPs and considering the timing of the transcriptional activation or inactivation. For this latter case, in order to have enough genes to compare for each case, the time breaks were grouped in order to represent embryonic (E15.5 Ngn3^eGFP/+), nascent (P1, P4), and older (P12, P60) β cells. The p-values resulting from the Fisher’s exact test represent the probability that the observed numbers of genes belonging to a specific case (e.g., “genes belonging to positive MPs” or “genes belonging to a positive MP and mainly activated in P60 β cells”) and annotated with a functional category/GSIS process have resulted from random sampling. Results showing a false discovery rate (FDR)–adjusted p-value lower than 5% were considered significant.

Antibody Staining, Ca²⁺ Imaging, in vitro GSIS, Immunoassays, and Transmission Electron Microscopy

Antibody staining followed routine protocols. Rabbit anti-amylase (Sigma–Aldrich) and guinea pig anti-insulin (Dako) were used at 1:1,000 dilution. Cy5-conjugated donkey anti–guinea pig and FITC-conjugated donkey anti-rabbit (Jackson Immunology) were used at 1:2,000 (Zhao et al., 2010). Hand-picked islets were used for whole-mount staining and imaged with confocal microscopy.

After 2-h recovery in RMPI1066 (5.6 mM glucose, 10% FBS), hand-picked islets were used for GSIS, immunoassays, and transmission electron microscopy (TEM). GSIS used basal glucose of 2.8 or 5.6 mM and stimulatory glucose of 20 mM as in Zhao et al. (2010). The percentage of insulin release from starting islets within a 45-min window was assayed. Immature and mature vesicles were classified and quantified according to electron density, with ImageJ. Ca²⁺ imaging followed protocols in Jacobson et al. (2010). All antibodies were from Jackson ImmunoResearch. For all the experimental comparisons, two-sided Student t-test was applied when the number of samples was >30; otherwise, Wilcoxon rank-sum test was used. For all the statistical tests, FDR-adjusted p-values < 5% were considered significant.

Results

Exploring Mip^eGFP Mice for β-Cell Isolation and Gene Expression Studies

Before applying the temporal analysis on both embryonic and postnatal stages, we first checked the utility of the Mip^eGFP cell–based gene expression for studying β-cell maturation at postnatal stages. Therefore, we examined whether eGFP expression labels all the β cells in the Mip^eGFP mice used for β-cell purification with FACS. As reported by Katsuta et al. (2012), the levels of eGFP in insulin⁺ cells greatly varied (Figures 1A,B), appearing as eGFP^High and eGFP^Low cells. Yet most of the insulin⁺ cells expressed detectable eGFP (1,143/1,216 counted insulin⁺ cells express eGFP; Figure 1A). These data suggested that collecting both eGFP^High and eGFP^Low cells (Figure 1B) will provide representative β cells in islets. We also examined whether the Mip^eGFP transgene interferes with endocrine islet function. At P4, P12, and P60, Mip^eGFP islets showed similar GSIS profiles as those of control islets (Figure 1C), suggesting the lack of detectable effects of Mip^eGFP transgene on GSIS.

FIGURE 1

Figure 1. Mip^eGFP islets maintain glucose responsiveness. Hand-picked, dissociated or intact, islets were used for these assays. (A) Immunostaining of P4 Mip^eGFP islet cells. Arrowheads: insulin⁺ cell with undetectable EGFP. Arrows: debris (judged by size) picking up secondary antibodies. (B) EGFP analysis in P4 Mip^eGFP islet cells via flow cytometry. R1 and R2 are the EGFP^Lo and EGFP^Hi cells collected for RNA-seq. (C) GSIS of islets from Mip^eGFP mice and wild-type littermates at P4, P12, and P60. Shown are the percentages of total insulin from starting islets released in 30 min. ^#p ≥ 0.1, Wilcoxon rank-sum test (n = 5, number of individual GSIS assays). At least six mice were used for each study, with islets from two to three mice mixed as technical and biological repeats.

As further validation of using Mip^eGFP cell–based gene expression for monitoring β-cell maturation, we examined the expression of 60 genes in Rip^mCherry β cells with RT-PCR and determined whether the results consistently match the RNA sequencing (RNA-seq) results in Mip^eGFPβ cells. Rip^mCherry mice express mCherry in β cells under the control of a rat Insulin 2 promoter and a SV40 polyA. These mice have normal β-cell function and GSIS (Zhu et al., 2015). We purified β cells to ∼98.4% purity from Rip^mCherry mice at P1, P4, P12, and P60 (Figure 2A). The set of 60 genes chosen for RT-PCR includes both genes known for β-cell maturation and genes that are not expressed in β cells as positive and negative controls, respectively. We also included genes required for proliferation, differentiation, β-cell electrical activity, vesicular biosynthesis and secretion, stress responses, and metabolism, because of their established roles in β-cell production and GSIS (Figure 2B and Supplementary Figure 1). We randomly picked candidates expressed at both high (such as Pax6, Pdx1, and Hsps that are known in β-cell differentiation and function) and low levels (e.g., Pax4, Ngn3, and Ptf1a, involved in progenitor differentiation and down-regulated in β cells) to better represent the data set. Hprt was used as an internal control for PCR.

FIGURE 2

Figure 2. Verification of RNAseq data with RT-PCR in P1-P60 beta cells. The gene expression in P1 beta cells was used as a reference point for other stages. (A) An example of sorted P4 beta cells from Rip^mCherry mice. Note the two mCherry-negative cells (arrowheads) in the sorted population. (B) The expression levels of 20 genes, assayed with RT-PCR in Rip^mCherry+ cells (black lines) and RNAseq in Mip^eGFP cells (blue lines). *p < 0.05, Wilcoxon rank-sum test. For each time point, three RNA preps [each from one (P60) to three mice (P1 to P12)] were used for RT assays.

The expressions of non–β-cell genes, such as Amy1 (Figure 2B), Arx, Gcg, and Ptf1a (Supplementary Figure 1), displayed disparate patterns between Mip^eGFP RNAseq and Rip^mCherry RT-PCR results. These findings are consistent with a possibility that exocrine cells cannot be removed at 100% efficiency, and different samples could result in unpredictable acinar contaminations. Indeed, the per-cell levels of Amy1 in P1 sorted β cells are between 0.5 and 3% of that in total pancreata, consistent with the high purity of the sorted β cells (with ct ∼21 in total pancreas vs. ct ∼27 in sorted β cells). Among the rest of the 55 genes, 33 showed no significant difference in expression at all stages (Figure 2B and Supplementary Figure 1); 17 genes showed similar trends of expression dynamics with one stage that displayed significant differences between the two data sets, including Hes1, Kcnj2, Kcnk9, Syt4, Ldha, Pfkl (Figure 2B), Dnajb9, Fboxo2, His2hc31, Irx2, Kcnj5, Kcnk3, Rab3a, Tmed5, Tpt1, Vgf, and Zcchc12 (Supplementary Figure 1); 5 genes, i.e., Enol (Figure 2B) and Atf3, Cbs, Hspa1b, and Syt14 (Supplementary Figure 1), displayed significant difference at two stages, disrupting the dynamic trends of gene expression (Supplementary Figure 1). Notably, all known genes involved in maturation (MafA, MafB, NeuroD, and Ucn3) showed identical expression dynamics between Mip^eGFP and Rip^mCherry β cells (Figure 2B). These combined findings suggest that most of the gene expression dynamics obtained from Mip^eGFP cells reflect that of wild-type β cells.

Temporal Transcriptome Analysis of β Cells and Progenitors From RNA-Sequencing Data

RNA-seq data were generated for six cell populations from endocrine progenitors (E15.5 Ngn3^eGFP/eGFP and E15.5 Ngn3^eGFP/+) to mature β cells monitoring the postnatal stages P1, P4, P12, and P60 (see section “Materials and Methods”). There were 37 million to 85 million raw reads produced per sample. An average of 84.5% of the reads were uniquely mapped. Normalized expression counts were obtained for 39,016 Ensembl genes. Genes with fewer than 10 counts in each stage on average across the biological replicates were filtered out, leaving 18,445 expressed genes. The expression values of these genes were then log2-transformed (zero counts were kept as 0) for further analyses.

The FunPat pipeline identified 4,682 differentially expressed genes (Supplementary Tables 1, 2) across the six cell populations. Forty-five of these genes, including proteases (amylase, trypsin, CPAs, etc.) and nucleases (RNAse 1), were highly expressed in acinar cells (unpublished data). We therefore suspect that these mRNA could come from acinar contamination when FACS can only achieve 98–99% cell purity for the β cells (Supplementary Figure 2A). Indeed, immunoassays showed that young islets could tightly associate with amylase⁺ cells, which make them impossible to remove by hand-picking. Moreover, we could not detect amylase proteins in β cells despite the substantial level of mRNA in sorted β cells (Supplementary Figures 2B,C and Supplementary Table 1). Therefore, these 45 genes were removed for further analysis.

Considering the sign of the differential expression over the baseline at their corresponding time breaks, the TPs identified from the remaining genes were summarized into 11 “negative” and 18 “positive” MPs. Overall, more genes related to negative (3,436) rather than to positive (1,201) MPs were identified, as early progenitors are more heterogeneous than maturing β cells. Figure 3 displays the resulting MPs grouped according to the most representative time break. Twenty-one genes clustered alone, and their temporal profiles are reported in Supplementary Figure 3.

FIGURE 3

Figure 3. Mature patterns (MPs) identified in maturing β cells. Positive and negative MPs classified according to their time break (highlighted in red) and the dynamic expression across all stages. Note that the line y = 0 represents the basal gene expression at E15.5 in Ngn3^eGFP/eGFP cells that do not express endocrine specific genes. The capital letters are referred to in the MP description in Supplementary Table 1 and in the main text.

Most genes stabilized their transcriptional activity between P1 and P12 (MPs “A,” “D,” “J,” “K,” “N,” “O,” “T,” “U,” and “C1” in Figure 3). Interestingly, 254 genes in MP “A” showed steady-state expression from E15.5 Ngn3^eGFP/+ to P60, suggesting that these genes are involved in both β-cell differentiation and function. Several dynamics associated with highest/lowest expression level from the baseline were observed at P60. The positive MPs “U-Y” include known genes regulating cellular stress (DNAj9 and Hsps), membrane permeability (Lrrc55), cargo transport (Ina and Neb), and signaling (Gdf3, Gria 2, Prlr, and RGS). On the other hand, the negative MPs “Z-C1” displayed various degrees of decreased expression between P12 and P60, including down-regulation of genes, like Slc16a1 and Ldha, shown to be necessary for β-cell maturation (Ainscow et al., 2000).

More dynamic changes were observed in some MPs with few genes. For example, the 25 genes in MPs “B, C” showed decreased expression after E15.5 Ngn3^eGFP/+. Their main positive roles are likely for endocrine differentiation, but not GSIS. Examples include Ghrl, Arx, and Irx2 (MP “C”). Ghrl suppresses GSIS, whereas Arx and Irx2 are determinants of α cells (Collombat et al., 2003; Petri et al., 2006). Negative MPs “G” and “I” showed decreased expression between Ngn3^eGFP/+ and P1, followed by a postnatal increase. These include Cbs, Gabra4, and Wdr86 that regulate metabolism and neurotransmission of vesicle fusion (MacDonald et al., 2005; Benitez et al., 2012), consistent with the importance of metabolism for GSIS. MP “P”, showing continuously increased expression between P1 and P60, includes Ucn3, shown to stimulate insulin secretion in functional β cells (Blum et al., 2012). MP “R”, showing increased expression between P1 and P12 followed by a decrease at P60, includes Sycn, a secretory granule protein acting as Ca²⁺-sensitive regulator of exocytosis (Li et al., 2007), and Reg1, involved in proliferation of β cells. Interestingly, Reg2 is among the 21 single-gene MPs, and it showed its highest expression at P12 as does Reg1 but with positive differential expression at Ngn3^eGFP/+ and P1 (Supplementary Figure 3). The functional implication of Reg1 in the maturation process is not clear yet.

Main Functional Terms Characterizing β-Cell Maturation

To better summarize the biological information, the functional terms selected by FunPat were grouped into 17 main functional categories (Table 1) related to common ancestor terms and a manually-defined list of reference GO ancestor terms (Supplementary Data Sheet 2). We then examined through Fisher’s exact test the enrichment of each category in both positive and negative MPs (Table 1) in endocrine progenitors, nascent (P1-P4), and older β cells (P12-P60). Overall, 8 of the 17 functional categories were significantly enriched in both the positive and negative MPs, including “Adhesion, communication, aggregation, and migration,” “binding,” “biological regulation and behavior,” “cellular component organization or biogenesis,” “environmental information processing, response to stimulus, and signaling,” “membrane,” “organelle,” and “unclassified.” These findings showed the co-presence of both positive and negative regulators for β-cell maturation. In the positive MPs, all the eight functional categories (1,558 genes) showed enrichment in P1-P4 β cells, and four of them (364 genes) showed further enrichment in mature P12-P60 β cells. These results suggest that most of the genes positively required for maturation reached their highest expression by P4. Only a small number of genes, 364, need to be expressed in later mature β cells. Interestingly, the negative MPs showed different enrichment profiles; all except the “unclassified” category showed enrichment in mature β cells (Table 1). These results suggest that the down-regulation of the genes preventing β-cell maturation is a slower process. Genes in these groups might represent potential limiting factors for maturation.

TABLE 1

Table 1. The 17 functional categories and their enrichment for either positive or negative main patterns (independent of time breaks).

Several functional categories were enriched in either the positive or negative MPs, respectively. “Localization and transport” and “system process” were enriched in the positive MPs. “Cell cycle, proliferation, growth, and death,” “biosynthesis and catalytical activity,” “localization and transport,” and “system process” were specifically enriched in negative MPs, with the most significant enrichments occurring in mature β cells. Identification of the “cell cycle, proliferation, growth, and death” category is consistent with the established studies showing the reduced proliferation in β cells.

Main GSIS Processes Characterizing β-Cell Maturation

Because of the broad definition of the GO-based functional annotation, the above analyses could not reveal some specific biochemical pathways directing β-cell maturation. We therefore further analyzed the data against manually defined 10 GSIS processes according to known pathways reported in the literature (MacDonald et al., 2005; Benitez et al., 2012). Results were reported in Table 2 and Supplementary Data Sheet 2. For the positive MPs, 6 of the 10 processes displayed significant enrichment in immature β cells, including “calcium-mediated processing,” “GTPase and G-protein activity,” “insulin processing and signaling,” “ion transport and homeostasis,” “membrane potential and ion channels,” and “vesicle-mediated transport and secretion by cell.” Only “calcium-mediated processing,” among other enriched processes, showed further enrichment in mature β cells. These data indicate that most of the genes positively needed for β-cell maturation reached their plateau by P4. Yet, their expression is not sufficient to define maturity, and “calcium-mediated processing” is likely the key molecular mechanism limiting the maturation process among the positive regulators.

TABLE 2

Table 2. Enrichment results on GSIS processes for positive and negative MPs, according to time breaks.

Three GSIS processes resulted enriched in the negative MPs. As expected from previous studies, “glycolysis, glucose processing, pyruvate metabolism, and TCA cycle” and “oxoacid metabolic process and fatty acid activity” belong to this group. Their reduction reached the lowest level by P4 (Table 2). Another process “protein tyrosine and serine/threonine kinase activity” was also enriched in the negative MPs, specifically in mature P12-P60 β cells. This finding is consistent with the notion that hormonal regulation is essential for β-cell function.

Protein–Protein Interaction Network of Calcium-Mediated Processing

We next focused on “calcium-mediated processing,” enriched for positive MPs with time breaks at both nascent and older β cells, form protein–protein interaction (PPI) networks describing potential limiting factors for β-cell maturation. Annotations available from STRING database (von Mering et al., 2003)⁴ were used. 109 of the 142 genes belonging to this process were able to form a PPI network with 383 interactions (Figure 4). 54 of the 109 interacting genes belong to positive MPs, mostly representing channels and G-proteins, consistent with their positive roles in GSIS. The remaining 55 genes assigned to negative MPs were found associated with growth factor signaling (transforming growth factor β and insulin-like growth factor), consistent with lack of significant proliferation of mature β cells. Only 13 genes displayed enriched expression at P12 and P60. These included some vesicular proteins known to be involved in exocytosis process and all linked by Vamp-2/synaptobrevin, whose expression level resulted mainly established at birth: syntaxin-1 (Syn1), the NMDA receptor Grin1, whose deletion was recently associated with a higher degree of islet GSIS (Marquard et al., 2015), and synaptotagmins Syt4 and Syt5 (Gauthier and Wollheim, 2008; Huang et al., 2018). Syt7 was also found associated with calcium-mediated processing, but it showed no PPIs and highest expression by P1 (Supplementary Table 1). Another complex of cadherins (Cdh4, Cdh7, and Cdh8), recently associated with the increase of GSIS activity in β cells (Parnaud et al., 2015), showed their highest transcriptional activity in older β cells.

FIGURE 4

Figure 4. PPI network of calcium-mediated processing. The 109 genes selected and associated with calcium-mediated processing by FunPat correspond to a PPI network of 383 interactions. Gray nodes: genes belonging to negative MPs. Yellow and red nodes: genes belonging to positive MPs. Red nodes: genes belonging to MPs with time breaks at P12 or P60.

Temporal Development of GSIS During β-Cell Maturation

To correlate gene expression dynamics with β-cell maturation, we examined the GSIS of islets at postnatal stages (Figures 5A,B). P1 islets showed higher basal insulin secretion compared to adult islets (2.8 and 5.6 mM glucose), whereas higher glucose did not enhance insulin secretion. P4, P12, and adult islets all significantly responded to high glucose (20 mM) with insulin secretion, but P4 islets resulted not mature, showing higher basal insulin secretion compared to P12/P60 islets. Interestingly, preincubating isolated islets at 5.6 mM glucose before GSIS assays eliminated insulin secretion at high glucose in P4 islets, but not in P12 or P60 islets. The presence of releasable insulin vesicles at low glucose before P12 is consistent with our transcriptional data, showing that most transcripts for insulin vesicle biosynthesis are established before maturation completes (Figure 5 and Supplementary Table 1).

FIGURE 5

Figure 5. Temporal development of GSIS during β-cell maturation. Hand-picked wild-type islets were used to recover for 2 h before GSIS and immunoassays. For GSIS, the percentages of total insulin from starting islets released in a 45-min window are shown. For all assays, n ≥ 4 independent GSIS assays were done. (A) Islet GSIS for islets preincubated 2.8 mM glucose for 1 h. (B) GSIS after islet preincubated in 5.6 mM glucose for 1 h. (C–F) Morphologies of typical isolated islets. Shown are single-confocal planes of whole-mount islets or pancreatic sections that were co-stained for insulin (red), glucagon (blue), and somatostatin (green). *FDR-adjusted p ≤ 0.05, Wilcoxon rank-sum test (n ≥ 5).

To examine whether GSIS of P1-P60 islets resulted from their intercellular communication with each other, we compared the morphology of hand-picked islets. For islet size, ∼30–100 islets per mice (P1 to P60) were scored. Photographs were taken under a stereoscope, measuring the diameters. At least six mice (three males/three females) were examined at each stage. For looking at the cell-type ratio, ∼20 microscopic fields were scored to examine the percentage of islet cells that express insulin at each stage. At least six mice were scored for each stage. For β cell size, the (insulin⁺ areas)/nucleus (DAPI) was scored with the same fields. The obtained results showed that the size and percentage of β cells increase as they age (Figures 5C–F). Yet they all displayed similar morphology, with β cells clustered in the central region and other cell types in the periphery. Consistently, the gene annotations linked to the functional category “Cell adhesion, communication, aggregation, and migration” are enriched for positive MPs reaching steady state within P4 (Table 2). Therefore, islet organization alone may not account for the different GSIS properties of islets at different ages. Neither does the islet cell-type composition, because P4 and P12 islets display similar β cell/endocrine ratios yet have different maturity.

Vesicle Biosynthesis in β Cells Correlates With the Temporal Expression of Vesicular Genes

Previous findings suggested that proper vesicular packaging contributes to β-cell maturation (Blum et al., 2012; Goodyer et al., 2012). As our data showed that transcripts for most vesicular components reached a plateau in their expression by P4 (Table 2 and Supplementary Table 1), we determined the vesicular density in β cells during their maturation process. Indeed, a few P1 β cells showed well-defined mature vesicles with dense-core insulin crystals (Figure 6A), whereas others had mostly immature vesicles with electron-light core (Figure 6B) consistent with the results previously reported in the literature (Blum et al., 2012). By P4, most β cells had mature vesicles of similar appearance (Figures 6C–E), in terms of both vesicle density and morphology (Figure 6F). These results, combined with the temporal sequence of GSIS development, suggest that producing morphologically normal vesicles precedes maturation by over a week. These findings are also consistent with the fact that post-transcriptional regulation does not prominently regulate the production of vesicular proteins, because the appearance of the vesicles structure temporally coincides with gene transcription up-regulation.

FIGURE 6

Figure 6. Normal-looking vesicles are present in immature β cells. Hand-picked islets were fixed and examined with TEM. (A,B) Two P1 β cells showing different status of vesicular packaging. (C–E) Representative images of P4, P12, and P60 β cells. Arrows point to normal-looking “mature” vesicles. Arrowheads point to “immature” vesicles with less dense cores. (F) The levels of mature and immature vesicles in β cells are established by P4 based on quantitation with at least 26 microscopic views (1–3 β cells in each view) counted. Images used were from different blocks of at least four mice at each stage. *p < 0.05, two-sided Student t-test (n > 30).

Glucose-Induced Ca²⁺ Influx Is Established in Immature β Cells

The PPI network of the calcium-mediated process also included several genes encoding channel proteins and metabolic enzymes reaching a plateau of expression before P4 (Figure 3 and Table 2). We therefore examined glucose-induced Ca²⁺ influx in islets at different stages. Because Ca²⁺ influx depends on proper glucose transport, metabolism, oxidative phosphorylation, and activation/inactivation of multiple channels, proper Ca²⁺ activity in β cells will likely reveal the production and assembly of all protein complexes and pathways involved in these processes.

P1 islets displayed very low glucose-induced Ca²⁺ influx (Figure 7A), although they showed recognizable oscillating patterns, a property of β cells (Figure 7B). P4 to P60 islets showed significantly higher glucose-induced Ca²⁺ influx than P1 islets, with β-cell specific oscillations (Figures 7C–E). These findings suggest that the molecular machineries for glucose transport, metabolism, ATP production, ATP-mediated blockage of K⁺ channels, and voltage-gated Ca²⁺ channels are already present in immature β cells at P4. As the mRNAs coding for the proteins involved in the above processes also reached a plateau at P4, these findings indicate that post-transcriptional regulation is not prominently involved in controlling the production of these proteins. Finally, a notable and recurrent difference between immature and mature islets was observed in the 1- to 2-min delay between the applications of high glucose to the Ca²⁺ influx in mature islets (P < 0.001, Figure 7A). The reason and significance of this delay are not clear yet.

FIGURE 7

Figure 7. Glucose-induced Ca²⁺ influx is established in immature β cells. Free Ca²⁺ was recorded by Fura2 fluorescence. Basal glucose was used for the first 200 s, and then 20 mM glucose was used to stimulate islets. (A) Average Ca²⁺ responses of islets at each stage (n ≥ 69). (B–E) Examples of oscillating Ca²⁺ in three different islets. The color of lines here indicate different islets of same stages. The numbers of mice used were as follows: P1: 10; P4: 10; P12: 10; P60: 4.

Discussion

Our multi-staged transcriptome analysis using FunPat pipeline revealed several candidate genes, dynamic trends, and biological processes able to distinguish genes required for β-cell differentiation (generating insulin⁺ cells) and/or maturation (insulin⁺ cells gaining proper glucose response). Compared with the several pairwise gene expression analyses (Jermendy et al., 2011; Blum et al., 2012; Dhawan et al., 2015), our studies confirmed the importance of proper metabolism, calcium signaling, and vesicular biosynthesis for β-cell maturation. The multi-stage dynamic analysis allowed us to reveal several previously unrecognizable features but also to confirm the expression patterns of several well-recognized markers for predicting β-cell maturation (Benitez et al., 2012) including MafA and Ucn3, found also in other recent studies at a single-cell level (Qiu et al., 2017; Augsornworawat et al., 2020).

Besides the nature of the transgene, we have considered several variables that could affect our gene expression, including the islet isolation, dissociation (Khan et al., 2016), and mechanical sorting (Beliakova-Bethell et al., 2014). All these processes can in theory alter the gene expression. Yet, our RNAseq-based results directly correlate with tissue-staining–based gene expression, such as MafA (Hang et al., 2014), NeuroD1 (Gu et al., 2010), and Ucn3 (Blum et al., 2012), suggesting that these technical issues will unlikely invalidate most of the gene expression patterns. It is worth considering that these findings are based on mouse models and that there are differences with human β cells. For example, MafB is expressed in human but not mouse β cells (Arda et al., 2016; Xin et al., 2016; Tritschler et al., 2017), whereas Ucn3 expression is much higher in mouse than human β cells (Xin et al., 2016). However, these observations have been made analyzing β cells by considering no more than two maturation stages; therefore, we believe that our study will be useful for future comparisons with temporal gene expression data from human β cells.

Our RNAseq-based studies could not address whether post-transcriptional regulation prominently regulates β-cell maturation. However, by examining the presence of the vesicular structures and biochemical pathways for insulin secretion and Ca²⁺ influx, we found that the appearance of required proteins closely matched the dynamic activation of gene transcription. These findings suggest that translational regulation is not a major regulatory mechanism for most of the genes involved in β-cell maturation and GSIS.

One of our main conclusions is that the expression of most maturation genes reaches a plateau before P4. Many of these gene products are involved in insulin biosynthesis, signal transduction of transmembrane receptors, vesicle transport, and calcium-mediated processing (Table 2, positive MPs). Indeed, producing Ca²⁺-related proteins or signal transducers is likely needed for proper stimulus-secretion coupling. Such transcription profiles correspond to detection of significantly higher density of mature vesicles (Figure 6F) and glucose-induced Ca²⁺ influx (Figure 7A) in P4 with respect to P1.

In addition, incubating P4 β cells with low basal glucose (2.8 mM, Figure 5A) led to significant insulin secretion when glucose stimulus was switched to 20 mM, but not when basal glucose levels were increased (5.6 mM, Figure 5B). This result suggests that repressing insulin secretion at low glucose could be a limiting step for β-cell maturation. Indeed, incubating immature β cells with higher basal glucose appears to deplete the releasable vesicle pool so that switching to higher glucose could no longer trigger further insulin secretion, as already reported (Blum et al., 2012). The implications of these observations are as follows: (1) vesicles in immature β cells, even if morphologically normal-looking (Figure 6), are not equally releasable, a feature that has been proven by real-time secretion assays (Hou et al., 2012; Hoboth et al., 2015); (2) insulin vesicles need to desensitize themselves from glucose-derived signals to abstain from releasing at basal glucose in the maturation process. This could theoretically be achieved by limiting the Ca²⁺ influx at basal glucose or modulate their Ca²⁺ sensitivity, or both. In support of this idea, enrichment analysis found only “calcium-mediated processing” to be significant among the GSIS processes for positive MPs characterizing both immature and mature β cells (Table 2). The corresponding MPs and PPI network (Figures 3, 4) highlight the interaction of synaptotagmins and cadherins at later stages of maturation, suggesting a possible role in controlling vesicle release and Ca²⁺ influx (Huang et al., 2018).

Finally, the transcriptome data analysis also highlighted the importance of negative GSIS regulation in β-cell maturation, generally represented by MPs showing an expression decrease between P12 and P60. Focusing on GSIS processes, enrichment analysis found only “protein tyrosine and serine/threonine kinases activity” to be significant for negative MPs at late maturation stages (Table 2), suggesting unappreciated phosphorylation processes linked with vesicle biosynthesis and release. This study can be considered a resource for further integrations/comparisons to human β-cell development from embryonic or induced-pluripotent stem cells (Weng et al., 2020). The complete database of TPs, genes, and functional terms available in Supplementary Tables 1, 2 will aid future studies to better characterize the regulatory role of these genes and critical steps for in vitro–derived functional β cells.

Data Availability Statement

The dataset presented in this study can be found at ArrayExpress repository (https://www.ebi.ac.uk/arrayexpress/) with the accession number E-MTAB-2266.

Ethics Statement

The animal study was reviewed and approved by the Vanderbilt University Animal Care and Use Committee.

Author Contributions

TS, GG, CS, and MM conceived the work and designed the experiments. TS, GG, and CS wrote the manuscripts, with help from all listed authors. TS, EM, BD, and CS performed all the bioinformatics analysis. CH, YX, and LP isolated the cells and performed the RNA preparations and islet characterization. PD and DJ performed the Ca²⁺ imaging. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by grants from NIDDK (DK065949, DK128710, and DK125696 for GG and UO1DK089523 for MM), JDRF (1-2009-371 for GG), PRAT 2010 (CPDA101217 to BC), and Fondazione Aldo Gini (Borsa Gini scholarship to TS). Flow Cytometry experiments were performed in the VUMC Flow Cytometry Shared Resource, supported by the Vanderbilt Ingram Cancer Center (P30 CA68485) and the Vanderbilt Digestive Disease Research Center (DK058404). Confocal and TEM imaging were performed with VUMC Cell Imaging Shared Resource (supported by NIH Grants CA68485, DK20593, DK58404, DK59637, and EY08126).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Chris Wright and Roland Stein for useful comments and discussion. TS thanks the programme “Dipartimenti di Eccellenza 2018–2022” (Project code D15D18000410001) specifically appointed to her affiliation, the Department of Medical Sciences of University of Torino, from the Italian Ministry for Education, University and Research (Ministero dell’Istruzione, dell’Università e della Ricerca – MIUR).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.648791/full#supplementary-material

Footnotes

1. ^ https://github.com/itmat/Normalization
2. ^ http://www.informatics.jax.org/function.shtml
3. ^ http://cpdb.molgen.mpg.de/MCPDB
4. ^ http://string-db.org/

References

Adam, P. A., Teramo, K., Raiha, N., Gitlin, D., and Schwartz, R. (1969). Human fetal insulin metabolismearly in gestation. Response to acutelevation of the fetal glucose concentration and placental tranfer of human insulin-I-131. Diabetes 18, 409–416. doi: 10.2337/diab.18.6.409

PubMed Abstract | CrossRef Full Text | Google Scholar

Ainscow, E. K., Zhao, C., and Rutter, G. A. (2000). Acute overexpression of lactate dehydrogenase-A perturbs beta-cell mitochondrial metabolism and insulin secretion. Diabetes 49, 1149–1155. doi: 10.2337/diabetes.49.7.1149

PubMed Abstract | CrossRef Full Text | Google Scholar

Arda, H. E., Li, L., Tsai, J., Torre, E. A., Rosli, Y., Peiris, H., et al. (2016). Age-dependent pancreatic gene regulation reveals mechanisms governing human beta cell function. Cell Metab. 23, 909–920. doi: 10.1016/j.cmet.2016.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene ontology: tool for the unification of biology. The gene ontology consortium. Nat. Genet. 25, 25–29. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

Augsornworawat, P., Maxwell, K. G., Velazco-Cruz, L., and Millman, J. R. (2020). Single-cell transcriptome profiling reveals β cell maturation in stem cell-derived islets after transplantation. Cell Rep. 32:108067. doi: 10.1016/j.celrep.2020.108067

PubMed Abstract | CrossRef Full Text | Google Scholar

Beliakova-Bethell, N., Massanella, M., White, C., Lada, S., Du, P., Vaida, F., et al. (2014). The effect of cell subset isolation method on gene expression in leukocytes. Cytometry A 5, 94–104. doi: 10.1002/cyto.a.22352

PubMed Abstract | CrossRef Full Text | Google Scholar

Benitez, C. M., Goodyer, W. R., and Kim, S. K. (2012). Deconstructing pancreas developmental biology. Cold Spring Harb. Perspect. Biol. 4:a012401. doi: 10.1101/cshperspect.a012401

PubMed Abstract | CrossRef Full Text | Google Scholar

Benitez, C. M., Qu, K., Sugiyama, T., Pauerstein, P. T., Liu, Y., Tsai, J., et al. (2014). An integrated cell purification and genomics strategy reveals multiple regulators of pancreas development. PLoS Genet. 10:e1004645. doi: 10.1371/journal.pgen.1004645

PubMed Abstract | CrossRef Full Text | Google Scholar

Benninger, R. K., and Piston, D. W. (2014). Cellular communication and heterogeneity in pancreatic islet insulin secretion dynamics. Trends Endocrinol. Metab. 25, 399–406. doi: 10.1016/j.tem.2014.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bliss, C. R., and Sharp, G. W. (1994). A critical period in the development of the insulin secretory response to glucose in fetal rat pancreas. Life Sci. 55, 423–427. doi: 10.1016/0024-3205(94)90053-1

CrossRef Full Text | Google Scholar

Blum, B., Hrvatin, S. S., Schuetz, C., Bonal, C., Rezania, A., and Melton, D. A. (2012). Functional beta-cell maturation is marked by an increased glucose threshold and by expression of urocortin 3. Nat. Biotechnol. 30, 261–264. doi: 10.1038/nbt.2141

PubMed Abstract | CrossRef Full Text | Google Scholar

Collombat, P., Mansouri, A., Hecksher-Sorensen, J., Serup, P., Krull, J., Gradwohl, G., et al. (2003). Opposing actions of Arx and Pax4 in endocrine pancreas development. Genes Dev. 17, 2591–2603. doi: 10.1101/gad.269003

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhawan, S., Tschen, S. I., Zeng, C., Guo, T., Hebrok, M., Matveyenko, A., et al. (2015). DNA methylation directs functional maturation of pancreatic beta cells. J. Clin. Invest. 125, 2851–2860. doi: 10.1172/JCI79956

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Camillo, B., Irving, B. A., Schimke, J., Sanavia, T., Toffolo, G., Cobelli, C., et al. (2012). Function-based discovery of significant transcriptional temporal patterns in insulin stimulated muscle cells. PLoS One 7:e32391. doi: 10.1371/journal.pone.0032391

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Camillo, B., Toffolo, G., Nair, S. K., Greenlund, L. J., and Cobelli, C. (2007). Significance analysis of microarray transcript levels in time series experiments. BMC Bioinform. 8(Suppl. 1):S10. doi: 10.1186/1471-2105-8-S1-S10

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Gauthier, B. R., and Wollheim, C. B. (2008). Synaptotagmins bind calcium to release insulin. Am. J. Physiol. Endocrinol. Metab. 295, E1279–E1286. doi: 10.1152/ajpendo.90568.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodyer, W. R., Gu, X., Liu, Y., Bottino, R., Crabtree, G. R., and Kim, S. K. (2012). Neonatal beta cell development in mice and humans is regulated by calcineurin/NFAT. Dev. Cell 23, 21–34. doi: 10.1016/j.devcel.2012.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, C., Stein, G. H., Pan, N., Goebbels, S., Hornberg, H., Nave, K. A., et al. (2010). Pancreatic beta cells require NeuroD to achieve and maintain functional maturity. Cell Metab. 11, 298–310. doi: 10.1016/j.cmet.2010.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, G., Dubauskaite, J., and Melton, D. A. (2002). Direct evidence for the pancreatic lineage: NGN3+ cells are islet progenitors and are distinct from duct progenitors. Development 129, 2447–2457.

Google Scholar

Gu, G., Wells, J. M., Dombkowski, D., Preffer, F., Aronow, B., and Melton, D. A. (2004). Global expression analysis of gene regulatory pathways during endocrine pancreatic development. Development 131, 165–179. doi: 10.1242/dev.00921

PubMed Abstract | CrossRef Full Text | Google Scholar

Hang, Y., Yamamoto, T., Benninger, R. K., Brissova, M., Guo, M., Bush, W., et al. (2014). The MafA transcription factor becomes essential to islet β-cells soon after birth. Diabetes 63, 1994–2005. doi: 10.2337/db13-1001

PubMed Abstract | CrossRef Full Text | Google Scholar

Hara, M., Wang, X., Kawamura, T., Bindokas, V. P., Dizon, R. F., Alcoser, S. Y., et al. (2003). Transgenic mice with green fluorescent protein-labeled pancreatic beta -cells. Am. J. Physiol. Endocrinol. Metab. 284, E177–E183. doi: 10.1152/ajpendo.00321.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Hellerstrom, C., and Swenne, I. (1991). Functional maturation and proliferation of fetal pancreatic beta-cells. Diabetes 40(Suppl. 2), 89–93. doi: 10.2337/diab.40.2.s89

PubMed Abstract | CrossRef Full Text | Google Scholar

Henquin, J. C. (2011). The dual control of insulin secretion by glucose involves triggering and amplifying pathways in beta-cells. Diabetes Res. Clin. Pract. 93(Suppl. 1), S27–S31. doi: 10.1016/S0168-8227(11)70010-9

CrossRef Full Text | Google Scholar

Hoboth, P., Muller, A., Ivanova, A., Mziaut, H., Dehghany, J., Sonmez, A., et al. (2015). Aged insulin granules display reduced microtubule-dependent mobility and are disposed within actin-positive multigranular bodies. Proc. Natl. Acad. Sci. U.S.A. 112, E667–E676. doi: 10.1073/pnas.1409542112

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, J. C., Min, L., and Pessin, J. E. (2009). Insulin granule biogenesis, trafficking and exocytosis. Vitam. Horm. 80, 473–506.

Google Scholar

Hou, N., Mogami, H., Kubota-Murata, C., Sun, M., Takeuchi, T., and Torii, S. (2012). Preferential release of newly synthesized insulin assessed by a multi-label reporter system using pancreatic beta-cell line MIN6. PLoS One 7:e47921. doi: 10.1371/journal.pone.0047921

PubMed Abstract | CrossRef Full Text | Google Scholar

Hrvatin, S., O’Donnell, C. W., Deng, F., Millman, J. R., Pagliuca, F. W., DiIorio, P., et al. (2014). Differentiated human stem cells resemble fetal, not adult, beta cells. Proc. Natl. Acad. Sci. U.S.A. 111, 3038–3043. doi: 10.1073/pnas.1400709111

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C., Walker, E. M., Dadi, P. K., Hu, R., Xu, Y., Zhang, W., et al. (2018). Synaptotagmin 4 regulates pancreatic β cell maturation by modulating the Ca2+ sensitivity of insulin secretion vesicles. Dev. Cell 45, 347–361.e5. doi: 10.1016/j.devcel.2018.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobson, D. A., Mendez, F., Thompson, M., Torres, J., Cochet, O., and Philipson, L. H. (2010). Calcium-activated and voltage-gated potassium channels of the pancreatic islet impart distinct and complementary roles during secretagogue induced electrical responses. J. Physiol. 588, 3525–3537. doi: 10.1113/jphysiol.2010.190207

PubMed Abstract | CrossRef Full Text | Google Scholar

Jermendy, A., Toschi, E., Aye, T., Koh, A., Aguayo-Mazzucato, C., Sharma, A., et al. (2011). Rat neonatal beta cells lack the specialised metabolic phenotype of mature beta cells. Diabetologia 54, 594–604. doi: 10.1007/s00125-010-2036-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamburov, A., Pentchev, K., Galicka, H., Wierling, C., Lehrach, H., and Herwig, R. (2011). ConsensusPathDB: toward a more complete picture of cell biology. Nucleic Acids Res. 39(Database issue), D712–D717. doi: 10.1093/nar/gkq1156

PubMed Abstract | CrossRef Full Text | Google Scholar

Katsuta, H., Aguayo-Mazzucato, C., Katsuta, R., Akashi, T., Hollister-Lock, J., Sharma, A. J., et al. (2012). Subpopulations of GFP-marked mouse pancreatic beta-cells differ in size, granularity, and insulin secretion. Endocrinology 153, 5180–5187. doi: 10.1210/en.2012-1257

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, D., Vasu, S., Moffett, R. C., Irwin, N., and Flatt, P. R. (2016). Islet distribution of Peptide YY and its regulatory role in primary mouse islets and immortalised rodent and human beta-cell function and survival. Mol. Cell Endocrinol. 436, 102–113. doi: 10.1016/j.mce.2016.07.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Kieffer, T. J. (2016). Closing in on mass production of mature human beta cells. Cell Stem Cell 18, 699–702. doi: 10.1016/j.stem.2016.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Konstantinova, I., Nikolova, G., Ohara-Imaizumi, M., Meda, P., Kucera, T., Zarbalis, K., et al. (2007). EphA-Ephrin-A-mediated beta cell communication regulates insulin secretion from pancreatic islets. Cell 129, 359–370. doi: 10.1016/j.cell.2007.02.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Kushner, J. A., MacDonald, P. E., and Atkinson, M. A. (2014). Stem cells to insulin secreting cells: two steps forward and now a time to pause? Cell Stem Cell 15, 535–536. doi: 10.1016/j.stem.2014.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, C. S., Perreault, N., Brestelli, J. E., and Kaestner, K. H. (2002). Neurogenin 3 is essential for the proper specification of gastric enteroendocrine cells and the maintenance of gastric epithelial cell identity. Genes Dev. 16, 1488–1497. doi: 10.1101/gad.985002

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., Chen, P., Vaughan, J., Lee, K. F., and Vale, W. (2007). Urocortin 3 regulates glucose- stimulated insulin secretion and energy homeostasis. Proc. Natl. Acad. Sci. U.S.A. 104, 4206–4211. doi: 10.1073/pnas.0611641104

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

MacDonald, P. E., Joseph, J. W., and Rorsman, P. (2005). Glucose-sensing mechanisms in pancreatic beta-cells. Philos. Trans. R Soc. Lond. B Biol. Sci. 360, 2211–2225. doi: 10.1098/rstb.2005.1762

PubMed Abstract | CrossRef Full Text | Google Scholar

Marquard, J., Otter, S., Welters, A., Stirban, A., Fischer, A., Eglinger, J., et al. (2015). Characterization of pancreatic NMDA receptors as possible drug targets for diabetes treatment. Nat. Med. 21, 363–372. doi: 10.1038/nm.3822

PubMed Abstract | CrossRef Full Text | Google Scholar

Nielsen, K., Kruhoffer, M., Orntoft, T., Sparre, T., Wang, H., Wollheim, C., et al. (2004). Gene expression profiles during beta cell maturation and after IL-1beta exposure reveal important roles of Pdx-1 and Nkx6.1 for IL-1beta sensitivity. Diabetologia 47, 2185–2199. doi: 10.1007/s00125-004-1578-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishimura, W., Kondo, T., Salameh, T., El Khattabi, I., Dodge, R., Bonner-Weir, S., et al. (2006). A switch from MafB to MafA expression accompanies differentiation to pancreatic beta-cells. Dev. Biol. 293, 526–539. doi: 10.1016/j.ydbio.2006.02.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Osundiji, M. A., and Evans, M. L. (2013). Brain control of insulin and glucagon secretion. Endocrinol. Metabol. Clin. North Am. 42, 1–14. doi: 10.1016/j.ecl.2012.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Pagliuca, F. W., Millman, J. R., Gurtler, M., Segel, M., Van Dervort, A., Ryu, J. H., et al. (2014). Generation of functional human pancreatic beta cells in vitro. Cell 159, 428–439. doi: 10.1016/j.cell.2014.09.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Parnaud, G., Lavallard, V., Bedat, B., Matthey-Doret, D., Morel, P., Berney, T., et al. (2015). Cadherin engagement improves insulin secretion of single human beta-cells. Diabetes 64, 887–896. doi: 10.2337/db14-0257

PubMed Abstract | CrossRef Full Text | Google Scholar

Petri, A., Ahnfelt-Ronne, J., Frederiksen, K. S., Edwards, D. G., Madsen, D., Serup, P., et al. (2006). The effect of neurogenin3 deficiency on pancreatic gene expression in embryonic mice. J. Mol. Endocrinol. 37, 301–316.

Google Scholar

Qiu, W. L., Zhang, Y. W., Feng, Y., Li, L. C., Yang, L., and Xu, C. R. (2017). Deciphering pancreatic islet β cell and α cell maturation pathways and characteristic features at the single-cell level. Cell Metab. 25, 1194–1205. doi: 10.1016/j.cmet.2017.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Rorsman, P., Arkhammar, P., Bokvist, K., Hellerstrom, C., Nilsson, T., Welsh, M., et al. (1989). Failure of glucose to elicit a normal secretory response in fetal pancreatic beta cells results from glucose insensitivity of the ATP-regulated K+ channels. Proc. Natl. Acad. Sci. U.S.A. 86, 4505–4509. doi: 10.1073/pnas.86.12.4505

PubMed Abstract | CrossRef Full Text | Google Scholar

Rozzo, A., Meneghel-Rozzo, T., Delakorda, S. L., Yang, S. B., and Rupnik, M. (2009). Exocytosis of insulin: in vivo maturation of mouse endocrine pancreas. Ann. N. Y. Acad. Sci. 1152, 53–62. doi: 10.1111/j.1749-6632.2008.04003.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanavia, T., Finotello, F., and Di Camillo, B. (2015). FunPat: function-based pattern analysis on RNA-seq time series data. BMC Genomics 16:S2. doi: 10.1186/1471-2164-16-S6-S2

PubMed Abstract | CrossRef Full Text | Google Scholar

Schiesser, J. V., and Wells, J. M. (2014). Generation of beta cells from human pluripotent stem cells: are we there yet? Ann. N. Y. Acad. Sci. 1311, 124–137. doi: 10.1111/nyas.12369

PubMed Abstract | CrossRef Full Text | Google Scholar

Shahjalal, H. M., Abdal Dayem, A., Lim, K. M., Tak-il, J., and Ssang-Goo, C. (2018). Generation of pancreatic β cells for treatment of diabetes: advances and challenges. Stem Cell Res. Ther. 9:355. doi: 10.1186/s13287-018-1099-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Tabar, V., and Studer, L. (2014). Pluripotent stem cells in regenerative medicine: challenges and recent progress. Nat. Rev. Genet. 15, 82–92. doi: 10.1038/nrg3563

PubMed Abstract | CrossRef Full Text | Google Scholar

Tritschler, S., Theis, F. J., Lickert, H., and Bottcher, A. (2017). Systematic single-cell analysis provides new insights into heterogeneity and plasticity of the pancreas. Mol. Metab. 6, 974–990. doi: 10.1016/j.molmet.2017.06.021

PubMed Abstract | CrossRef Full Text | Google Scholar

von Mering, C., Huynen, M., Jaeggi, D., Schmidt, S., Bork, P., and Snel, B. (2003). STRING: a database of predicted functional associations between proteins. Nucleic Acids Res. 31, 258–261. doi: 10.1093/nar/gkg034

PubMed Abstract | CrossRef Full Text | Google Scholar

Weng, C., Xi, J., Li, H., Cui, J., Gu, A., Lai, S., et al. (2020). Single-cell lineage analysis reveals extensive multimodal transcriptional control during directed beta-cell differentiation. Nat. Metab. 2, 1443–1458. doi: 10.1038/s42255-020-00314-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Xin, Y., Kim, J., Okamoto, H., Ni, M., Wei, Y., Adler, C., et al. (2016). RNA sequencing of single human islet cells reveals type 2 diabetes genes. Cell Metab. 24, 608–615. doi: 10.1016/j.cmet.2016.08.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, A., Ohara-Imaizumi, M., Brissova, M., Benninger, R. K., Xu, Y., Hao, Y., et al. (2010). Galphao represses insulin secretion by reducing vesicular docking in pancreatic beta-cells. Diabetes 59, 2522–2529. doi: 10.2337/db09-1719

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, X., Hu, R., Brissova, M., Stein, R. W., Powers, A. C., Gu, G., et al. (2015). Microtubules negatively regulate insulin secretion in pancreatic beta cells. Dev. Cell 34, 656–668. doi: 10.1016/j.devcel.2015.08.020

PubMed Abstract | CrossRef Full Text | Google Scholar