The Genetic Program of Pancreatic β-Cell Replication In Vivo.

The molecular program underlying infrequent replication of pancreatic β-cells remains largely inaccessible. Using transgenic mice expressing green fluorescent protein in cycling cells, we sorted live, replicating β-cells and determined their transcriptome. Replicating β-cells upregulate hundreds of proliferation-related genes, along with many novel putative cell cycle components. Strikingly, genes involved in β-cell functions, namely, glucose sensing and insulin secretion, were repressed. Further studies using single-molecule RNA in situ hybridization revealed that in fact, replicating β-cells double the amount of RNA for most genes, but this upregulation excludes genes involved in β-cell function. These data suggest that the quiescence-proliferation transition involves global amplification of gene expression, except for a subset of tissue-specific genes, which are "left behind" and whose relative mRNA amount decreases. Our work provides a unique resource for the study of replicating β-cells in vivo.

Introduction

Pancreatic β-cells residing in the islets of Langerhans are highly specialized for secreting insulin in response to small elevations of blood glucose, thereby maintaining systemic glucose homeostasis. Reduced β-cell mass is a central feature of type 1 and type 2 diabetes, and strategies to enhance β-cell mass are sought as potential regenerative therapies for diabetes (1,2). Despite being a classic terminally differentiated cell type, β-cells are not entirely postmitotic, and rare β-cell divisions are responsible for homeostatic maintenance of β-cell mass (3–5). Furthermore, studies in rodents have shown that increased demand for insulin due to reduced insulin sensitivity or reduced β-cell mass triggers compensatory β-cell replication, leading to increased β-cell mass (6–8). The pathways regulating β-cell replication have been intensely investigated. Glucose acting via glucose metabolism and the unfolded protein response has emerged as a key driver of cell cycle entry of quiescent β-cells (7,9–11), and other growth factors and hormones have also been implicated in basal and compensatory β-cell replication under distinct conditions such as pregnancy (12,13). Intracellularly, multiple signaling molecules are involved in the transmission of mitogenic stimuli to the cell division cycle machinery of β-cells (14,15), including most notably the calcium–calcineurin–nuclear factor of activated T cells (NFAT) pathway (16,17). Despite this important progress, major gaps remain in our understanding of β-cell replication. To a large extent, this results from the rarity of β-cell replication in vivo. Dividing β-cells have been identified, counted, and localized in situ using immunostaining for markers such as Ki67 or BrdU, but it has not been possible so far to systematically characterize their biology on a genome-wide scale. Consequently, the molecular changes occurring in replicating β-cells in vivo remain largely unknown.

We have reasoned that the study of the transcriptome of replicating β-cells will require the isolation of these cells from islets while still alive so as to preserve their labile mRNA. To achieve this, we previously developed a transgenic mouse strain in which all cells express a fusion between the destruction box of cyclin B1 and green fluorescent protein (CcnB1-GFP) (18,19). In these mice, GFP is constitutively transcribed in all tissues from a weak phosphoglycerate kinase (PGK) promoter, but the GFP is rapidly degraded by the proteasome in quiescent cells. Degradation of the CcnB1-GFP fusion ceases when cells pass through the S-G2-M phases of cell cycle, thus increasing the level of GFP and enabling the identification and sorting of live replicating cells. We initially used this strain to characterize the transcriptome of replicating hepatocytes, providing a proof of concept as well as insights into the biology of cell division in the liver (19). In the current study, we used this transgene to isolate replicating and quiescent pancreatic β-cells from juvenile mice and analyzed their transcriptome by RNA sequencing (RNA-seq). Our data provide the first comprehensive view of the genes that are activated in spontaneously replicating β-cells in vivo, including many that have not previously been associated with cell cycle regulation in general and in β-cell proliferation in particular. In addition, we uncover a differential regulation of genes involved in β-cell function, which escape the global transcriptional increase in replicating cells.

Research Design and Methods

Mice

CcnB1-GFP mice, generated by infecting single-cell embryos with a lentivirus constitutively expressing a fusion between the destruction box of cyclin B1 and GFP (18), have previously been described (19). These mice were bred with RIP-Cre (20) and Rosa-LSL-Tomato (obtained from The Jackson Laboratory) mice. All animal experiments were performed in accordance with guidelines established by the joint ethics committee (institutional animal care and use committee) of The Hebrew University and Hadassah Medical Center. The Hebrew University is an Association for Assessment and Accreditation of Laboratory Animal Care International–accredited institute.

Flow Cytometry and Cell Sorting

Islets were isolated after collagenase perfusion as previously described (21). Islet cells were dissociated to single cells by treatment with trypsin/EDTA, stained with the Zn2+-binding dye N-(6-Methoxy-8-quinolyl)-p-toluenesulfonamide (TSQ) (22), and sorted on a FACS Aria (Becton Dickinson). For flow cytometry analysis, islet cells were fixed and permeabilized using Cytofix/Cytoperm and Perm/Wash solutions (BD Pharmingen). Immunostaining was performed in Perm/Wash solution using the following primary antibodies: guinea pig anti-insulin (Dako), goat anti-Pdx1 (a generous gift from Chris Wright [Vanderbilt University Medical Center]), and mouse anti-Nkx6.1 (Beta Cell Biology Consortium). Secondary antibodies were from Jackson ImmunoResearch Laboratories. Immunostaining for Ki67 was performed with the Alexa Fluor 647–conjugated mouse anti-Ki67 (BD Pharmingen).

RNA-seq and Analysis

Total RNA from ∼1,000–10,000 sorted β-cells was isolated by TRIzol (Invitrogen) extraction followed by an RNeasy Plus Micro kit (Qiagen). Poly(A)+-RNA was amplified using a MessageAmp II aRNA Amplification kit (Ambion) as previously described (23). Amplified RNA was fragmented and processed for library generation using the Illumina TruSeq RNA-Seq Sample Prep kit following the manufacturer's instructions. The RNA-seq libraries were subsequently sequenced on Illumina HiSEq. 2500. The resulting data were loaded into Galaxy (24). TopHat 2.0.5 was used to map the 50–base pair single reads to the NCBI37/mm9 assembly of the mouse genome. Gene expression analysis of the RNA-seq data was conducted using the R package DESeq (25). Differentially expressed genes between replicating and nonreplicating cells were selected based on false discovery rate (FDR) (Benjamini-Hochberg) ≤0.05. Gene set enrichment analysis based on the hypergeometric test was performed using the Genomica software (http://genomica.weizmann.ac.il/). In our analysis, we considered gene sets with a P value <0.01 and an FDR <0.05 to be significant. Gene sets were compiled from the MSigDB collections (26) or from the literature as described in the text.

Quantitative Real-Time PCR

Total RNA (1 ng) was used for first-strand cDNA synthesis using random primers (Roche, Indianapolis, IN) and reverse transcriptase (ImProm-II; Promega, Madison, WI). Real-time PCR was performed with SYBR Green Fast mix (Quanta Biosciences) in 96-well plates using the Bio-Rad real-time thermal cycler CFX96. All reactions were performed in triplicates with three biological replicates. The relative amount of mRNA was calculated using the comparative CT method after normalization to β-actin. Briefly, we calculated ΔCt values between each gene and β-actin, and ΔΔCt values were calculated between the ΔCt of each replicate and the average ΔCt for GFP-negative replicates. We used the following primers: Actb forward (5′-cacagcttctttgcagctcct-3′) and reverse (5′-gtcatccatggcgaactgg-3′), Top2A forward (5′-agcagattagcttcgtcaacagc-3′) and reverse (5′-acatgtctgccgcccttaga-3′), Ki67 forward (5′-ttgaccgctcctttaggtatgaa-3′) and reverse (5′-ttccaagggactttcctgga-3′), Pnlip forward (5′- ccacattgctggggaggctgg-3′) and reverse (5′-tcagcggggtccaaccctgt-3′), Hars forward (5′-ggagtggagcggatcttttc-3′) and reverse (5′-tcagtctctcctccagcagc-3′), Pabpc1 forward (5′-atcctctccatccgggtct-3′) and reverse (5′-ggtgtccaaagcacgttcc-3′), Slc38a4 forward (5′-aagcccagggtgctgagaa-3′) and reverse (5′-ggccgtccgggaattg-3′), Glul forward (5′-cccaacaagctggtgctatg-3′) and reverse (5′-ggttgctcaccatgtccatt-3′), Ero1lb forward (5′-tgattcgcaggaccactttt-3′) and reverse (5′-cccttgtagccagtgtaccg-3′), Slc30a8 forward (5′-tttgtggttgtcatcgaggc-3′) and reverse (5′-gtcaccacccagatgcaaag-3′), Gck forward (5′-cacggtggccacaatgatc-3′) and reverse (5′-cagccggtgcccacaa-3′), Cpe forward (5′-gtcacctcggctaaggatgg-3′) and reverse (5′-gtttccaggagcaagcaatcg-3′), Scg5 forward (5′-ggaagcggaggagtgtcaat-3′) and reverse (5′-tgccacaacattgtccaacc-3′), and Slc2a2 forward (5′-aaccgggatgattggcatgt-3′) and reverse (5′-ggcgaatttatccagcagca-3′).

Single-Molecule RNA FISH

Dissociated islet cells from Ccnb1-GFP mice were seeded on collagen-coated coverslips, fixed with 4% paraformaldehyde, and permeabilized with 70% ice-cold ethanol. We performed single-molecule RNA in situ hybridization (smRNA-FISH) as previously described (27). We used GFP fluorescence or smFISH for Ki67 mRNA (28) to specifically label replicating cells. Supplementary Table 5 lists probe sequences. Images were taken with a Nikon Ti-E inverted fluorescence microscope equipped with a ×100 oil immersion objective and a Photometrics Pixis 1024 CCD camera using MetaMorph software (Molecular Devices, Downingtown, PA). The image-plane pixel dimension was 0.13 μm. Cells were manually segmented, and dots were automatically counted using the ImageM software (29). For each cell, dots were counted only in the Z-stacks for which the cell was in focus. In addition, expression per cell was computed as the number of mRNA dots divided by the number of Z-stacks quantified per cell to include cells that only partially appeared within the image stack.

Mitochondrial Membrane Potential

Islets isolated from Ccnb1-eGFP mice were incubated overnight in standard RPMI 1640 culture medium (Gibco). The islets were then dissociated to single cells using Trypsin/EDTA. Single-cell suspensions were preincubated in Krebs-Ringer buffer supplemented with 0.5% BSA and 3 mmol/L glucose at 37°C for 30 min. After preincubation, cells were incubated for 90 min at 37°C in Krebs-Ringer buffer with either 3 mmol/L or 16 mmol/L glucose and tetramethylrhodamineethylester (TMRE) (7 nmol/L) (Life Technologies). Cell suspensions were analyzed using FACS Aria (Becton Dickinson).

Results

Sorting Live Replicating β-Cells

To define the transcriptional program of replicating β-cells, we used CcnB1-GFP mice to sort replicating and quiescent cells isolated from the islets of Langerhans of 1-month-old mice. A BrdU pulse-labeling experiment in vivo confirmed that the rare replicating β-cells were GFP positive (Supplementary Fig. 1). Sorting replicating β-cells proved particularly challenging given the low percentage of proliferating β-cells at this stage (∼1% GFP+ cells). In addition, in preliminary experiments we found that sorted GFP+ islet cells had large amounts of acinar-specific genes, reflecting contamination by replicating acinar cells and making it difficult to identify β-cell–specific changes in gene expression (not shown). To overcome this problem we used two different strategies to isolate β-cells prior to the split to GFP+ and GFP− populations. In one approach, we stained dispersed CcnB1-GFP islet cells with the β-cell–specific Zinc-binding fluorescent dye TSQ (22) and sorted TSQ+ cells (highly enriched for β-cells) that were GFP+ or GFP− (Supplementary Fig. 1). In an alternative approach, we generated compound insulin-Cre; Rosa-LSL-Tomato; CcnB1-GFP transgenic mice (Supplementary Fig. 1). From these mice we sorted Tomato+ cells (highly enriched for β-cells) that were GFP+ or GFP−. Each experiment included pooled islets from six to eight mice and was performed in biological duplicates, resulting in two sets of Tomato+ and two sets of TSQ+ sorted cell fractions, each set consisting of quiescent (GFP−) and replicating (GFP+) cells. Thus, altogether we obtained four independent samples of GFP+ cells and four matched samples of GFP− cells. To assess the purity of RNA preparations, we measured by real-time RT-PCR the RNA levels of the acinar-specific gene pancreatic lipase (Pnlip). Both TSQ and the Tomato reporter dramatically reduced the presence of contaminating acinar mRNA compared with unstained cells from dissociated islets (Fig. 1), suggesting a strong enrichment for β-cells. Real-time RT-PCR revealed an >80-fold increase in the mRNA level of the cell cycle marker Ki67 in the GFP+ cell fractions compared with unstained or GFP− cells. Together these results indicated that the two sorting strategies successfully enriched for replicating β-cells.

Figure 1

Ccnb1-GFP marks replicating cells and allows for FACS sorting of quiescent and replicating β-cells. A: qRT-PCR analysis of sorted unstained, TSQ-stained, and Tomato+ islet cells. Results were normalized to Actb. Pnlip and Ki67 expression levels are presented relative to the expression in unstained islets and GFP− sorted cells, respectively. B: GO enrichment analysis of upregulated genes (FDR ≤0.05) in replicating β-cells. Overrepresentation was computed using a statistical test based on hypergeometric distribution. Gene sets with a P value <0.01 and an FDR <0.05 were considered significant.

Genes Upregulated in Replicating β-Cells

The GFP+ fractions consisted of ∼1,000 cells each and yielded ∼5 ng mRNA, necessitating mRNA amplification prior to sequencing. We amplified RNA using T7-based in vitro transcription (23), prepared libraries, and performed RNA-seq. After evaluation of sequence quality and mapping to the mouse genome we used DESeq (25) to measure differential expression between replicating and quiescent β-cells (four biological replicates in each group). A heat map of sample-to-sample distance matrix confirmed that the TSQ and Tomato samples were highly similar, whereas gene expression differences between the GFP− and GFP+ samples were greater (Supplementary Fig. 1). This analysis revealed 3,312 differentially expressed genes (FDR ≤0.05), including 2,067 genes that were upregulated and 1,245 genes that were downregulated in replicating β-cells. The levels of the majority of transcripts detected in β-cells did not change between GFP− and GFP+ cells (but this depended on the method of normalization [see below]).

The genes most significantly induced in replicating β-cells included known cell cycle–regulated genes, in particular, genes involved in mitotic processes: cyclins, cyclin-dependent kinases, and genes encoding components of the cytoskeleton, microtubules, centrosomes, and mitotic spindle, as expected from cells in the S-G2-M phases of the cell cycle (Table 1). Gene ontology (GO) analysis of the proliferating β-cell signature revealed the expected enrichment for genes involved in DNA synthesis and mitosis and in proliferation-associated functions including DNA repair and chromatin remodeling (Fig. 1). We also observed upregulation of gene sets associated with RNA processing and splicing (Fig. 1 and Supplementary Table 1). These gene sets were described previously in the context of hepatocyte proliferation (19), but their link to the cell division cycle remains poorly understood. Importantly, some of the activated genes have not been previously annotated as belonging to the generic proliferation-associated GO categories and were not identified as cell cycle regulated in genome-wide expression studies (Table 2) and, hence, represent novel cell cycle–related genes.

Table 1

Top 30 genes upregulated in replicating β-cells, sorted from lowest FDR (Benjamini-Hochberg)

Gene symbol	Gene description	log2(FC)	FDR
Top2a	Topoisomerase (DNA) II α	7.4	3.254E-137
Cdk1	Cyclin-dependent kinase 1	7.4	5.773E-135
Hmmr	Hyaluronan-mediated motility receptor (RHAMM)	7.3	1.779E-125
Pbk	PDZ-binding kinase	6.9	3.012E-125
Cenpf	Centromere protein F	7.0	1.478E-123
Iqgap3	IQ motif containing GTPase-activating protein 3	8.3	1.26E-122
Nusap1	Nucleolar and spindle–associated protein 1	7.2	2.723E-122
Arhgef39	ρ Guanine nucleotide exchange factor (GEF) 39	8.2	1.213E-120
Cdca3	Cell division cycle–associated 3	6.8	4.879E-118
Cdca8	Cell division cycle–associated 8	7.0	4.879E-118
Neil3	Nei-like 3 (Escherichia coli)	7.6	6.64E-118
Cdc20	Cell division cycle 20	6.9	4.026E-117
Aurkb	Aurora kinase B	7.7	1.448E-116
Plk1	Polo-like kinase 1	7.7	1.447E-115
Tpx2	TPX2, microtubule-associated protein homolog (Xenopus laevis)	6.7	5.782E-115
Troap	Trophinin-associated protein	8.0	2.113E-114
Fam64a	Family with sequence similarity 64, member A	7.3	8.469E-114
Ckap2l	Cytoskeleton-associated protein 2–like	6.7	1.681E-113
Cenpa	Centromere protein A	6.6	2.528E-113
Cenpe	Centromere protein E	6.9	5.374E-113
Birc5	Baculoviral IAP repeat-containing 5	6.7	2.244E-112
Ccnb2	Cyclin B2	6.5	9.707E-112
Ska1	Spindle and kinetochore–associated complex subunit 1	8.2	9.707E-112
Ccna2	Cyclin A2	6.5	1.469E-111
Mxd3	Max dimerization protein 3	7.2	2.387E-111
Casc5	Cancer susceptibility candidate 5	7.5	5.009E-108
Aspm	Abnormal spindle–like, microcephaly associated (Drosophila)	6.7	8.609E-106
Ankle1	Ankyrin repeat and LEM domain–containing 1	7.5	9.805E-106
Kif23	Kinesin family member 23	6.8	9.805E-106
Shcbp1	Shc SH2-domain–binding protein 1	6.9	8.241E-104

DESeq was used to normalize the counts and compute log twofold change [log2(FC)] and FDR.

Table 2

Top 30 most significantly induced genes in replicating β-cells not previously annotated as cell cycle regulated and not present in cell cycle–related gene sets

Gene name	Gene description	log2(FC)	FDR
Sapcd2	Suppressor APC domain–containing 2	6.8	8.4E-102
D17H6S56E-5	DNA segment, Chr17, human D6S56E5	5.4	2.37E-91
Tcf19	Transcription factor 19	5.2	6.63E-81
4930579G24Rik	RIKEN cDNA 4930579G24	5.6	2.97E-78
BC030867	cDNA sequence BC030867	7.2	6.79E-70
Slc9a5	Solute carrier family 9 (sodium/hydrogen exchanger), member 5	6.3	1.24E-59
Mtfr2	Mitochondrial fission regulator 2	5.9	1.61E-59
Sez6	Seizure-related gene 6	6.9	2.4E-59
Trim59	Tripartite motif–containing 59	4.3	2.11E-57
Pbx4	Pre–B-cell leukemia homeobox 4	5.2	2.74E-55
Depdc1a	DEP domain–containing 1a	6.6	2.62E-54
4930427A07Rik	RIKEN cDNA 4930427A07 gene	4.2	8.27E-53
Ccdc18	Coiled-coil domain–containing 18	6.7	8.68E-52
Ccdc34	Coiled-coil domain–containing 34	3.9	1.71E-50
Cyp39a1	Cytochrome P450, family 39, subfamily a, polypeptide 1	4.0	4.44E-48
Rgs12	Regulator of G-protein signaling 12	4.6	5.37E-47
Zfp367	Zinc finger protein 367	3.5	3.44E-43
4930452B06Rik	RIKEN cDNA 4930452B06 gene	5.2	3.43E-42
2700094K13Rik	RIKEN cDNA 2700094K13 gene	3.3	5.08E-41
Zgrf1	Zinc finger, GRF-type–containing 1	3.9	1.46E-39
Katnal2	Katanin p60 subunit A-like 2	3.9	2.51E-37
Serpinb8	Serine (or cysteine) peptidase inhibitor, clade B, member 8	4.8	1.1E-35
Tbc1d31	WD repeat domain 67	3.3	1.94E-35
BC055324	cDNA sequence BC055324	4.5	2.21E-33
Ccdc77	Coiled-coil domain–containing 77	3.0	7.38E-30
Rgs3	Regulator of G-protein signaling 3	3.2	3.04E-28
Heatr7b1	Maestro heat-like repeat family member 2A	3.5	8.15E-27
Dnph1	2'-deoxynucleoside 5′-phosphate N-hydrolase 1	4.2	3.6E-26
Lrdd	Leucine-rich and death domain–containing	3.3	4.54E-26
Neurl1b	Neuralized homolog lb (Drosophila)	5.3	5.67E-26

List of 30 genes compiled from cell cycle–GO categories, Whitfield cell cycle, and Benporath cycling genes from MSigDB. Genes associated with cancer (candidate oncogenes, tumor suppressor genes, or genes highly expressed in cancer) appear in boldface type. DESeq was used to normalize counts, compute log twofold change [log2(FC)], and calculate FDR (Benjamini-Hochberg).

Downregulation of Genes Encoding β-Cell Function in Replicating β-Cells

We next examined the genes that were repressed in replicating β-cells. GO analysis of downregulated gene sets in replicating β-cells showed a striking association with biological processes critical for β-cell function, such as golgi vesicle transport, secretion, and regulation of insulin secretion (Fig. 2). Among these were genes involved in insulin processing (Ero1-like beta [Saccharomyces cerevisiae] [Ero1lb], proprotein convertase subtilisin/kexin type 2 [Pcsk2], and carboxypeptidase E [Cpe]) and secretion (secretogranin V [Scg5], vesicle-associated membrane protein 4 [Vamp4], syntaxin 4a [Stx4a], syntaxin 5a [Stx5a], and synaptotagmin-like 4 [Sytl4]), and genes encoding insulin granule proteins (RAB3A, member RAS oncogene family [Rab3a]; heat shock 70 kDa protein 5 [glucose-regulated protein, 78 kDa] [Hspa5]; cathepsin D [Ctsd]; and clusterin [Clu]) (30) (Table 3). Genes involved in glucose uptake and metabolism, including some that play a role in glucose-stimulated insulin secretion (GSIS), were also repressed [solute farrier family 2 (facilitated glucose transporter), member 2 (Slc2a2); glucose-6-phosphatase, catalytic, 2 (G6pc2); isocitrate dehydrogenase 1 (NADP+), soluble (Idh1); solute carrier family 25 (mitochondrial carrier; citrate transporter), member 1 (Slc25a1); and glutamate dehydrogenase 1 (Glud1). Interestingly, urocortin 3 (Ucn3), a marker of β-cell maturation (31), was also significantly downregulated in replicating β-cells. These results suggest that β-cell replication is associated with reduced expression of genes that are required for the main function of differentiated β-cells: production, processing, storage, and release of insulin in response to glucose. This is consistent with our observation that replicating hepatocytes transiently revert to a less mature, or differentiated, transcriptome (19). Repression of genes encoding β-cell function occurred also in more rapidly proliferating and less mature β-cells, isolated from mice at postnatal day 7 (Supplementary Fig. 4).

Figure 2

Genes repressed in replicating β-cells are involved in β-cell function and maturation. A: GO analysis of genes downregulated in replicating β-cells. B: Significant overrepresentation of genes upregulated in development in genes repressed in replicating β-cells (FDR ≤0.05) and, conversely, of genes repressed in development in genes induced in replicating β-cells (FDR ≤0.05). C: Genes induced in maturation and repressed in replication are associated with β-cell function. Genes repressed in maturation and activated in replication are associated with cell cycle biological processes. Overrepresentation was computed using a statistical test based on hypergeometric distribution.

Table 3

Genes involved in secretion, encoding for insulin granule proteins, and regulating GSIS and repressed in replicating β-cells

Secretion genes			Insulin granule genes			Genes regulating GSIS
Gene	log2(FC)	FDR	Gene	log2(FC)	FDR	Gene	log2(FC)	FDR
Sec23b	−0.73	7.80E-03	Rab3a	−0.83	1.38E-02	Slc2a2	−1.38	1.25E-08
Uso1	−1.19	1.4E-06	Ctsd	−0.61	4.15E-02	G6pc2	−1.83	7.85E-15
Trappc4	−0.67	3.01E-02	Cpn1	−1.26	9.54E-03	Idh1	−1.92	1.50E-13
Scamp1	−1.10	1.78E-05	Ctsl	−1.65	1.24E-11	Slc25a1	−0.69	2.10E-02
Osbpl5	−1.05	1.40E-03	Scg5	−1.06	8.27E-03	Glud1	−0.86	2.11E-03
Syn2	−0.90	1.65E-03	Hspa5	−0.61	3.05E-02	Ucn3	−0.67	2.23E-02
Rab3a	−0.83	1.38E-02	Tcn2	−1.24	2.71E-06
Mcfd2	−0.94	9.96E-03	Pcsk2	−0.72	9.36E-03
Tpd52	−0.70	1.34E-02	Sytl4	−0.76	5.06E-03
Rab3d	−0.93	7.59E-03	Pdia3	−0.66	1.86E-02
Sytl4	−0.76	5.06E-03	Clu	−0.74	1.16E-02
Arfgap3	−0.74	8.03E-03	Cpe	−0.82	1.99E-03
Stx4a	−0.69	3.40E-02	Aldoa	−2.83	1.46E-31
Stxbp3a	−1.10	5.88E-05
Blzf1	−1.01	6.95F-04
Srp54c	−0.94	2.62E-03
Sec23a	−0.94	1.62E-03
Sar1b	−0.83	2.66E-03
Snap23	−1.12	3.80E-05
Copb2	−0.95	2.40E-04
Scg5	−1.06	8.27E-03
Unc13b	−0.64	3.77E-02
Scfd1	−1.09	1.65E-05
Sar1a	−0.77	4.18E-03
Gars	−0.61	3.21E-02
Rab1	−0.65	2.34E-02
Cog1	−0.63	3.09E-02
Rab27a	−0.96	2.08E-04
Myo6	−0.59	4.56E-02

FDR (Benjamini-Hochberg) and log twofold change [log2(FC)] as calculated by DESeq after count normalization.

A recent study has revealed that the unfolded protein response (UPR) during mild endoplasmic reticulum stress correlates with, and triggers β-cell proliferation in, the context of high glucose (11). Strikingly, the UPR markers Hspa5/BiP and activating transcription factor 6 (Atf6) were repressed in replicating β-cells (Table 3 and data not shown). This finding suggests that although high glucose–mediated mild ER stress can stimulate cell cycle entry of β-cells, active UPR is not an integral part of the replication program of these cells at homeostasis.

To better define the transcriptome changes in replicating β-cells and delineate the relationship between β-cell replication and maturation, we generated gene sets from published transcriptomes of β-cells at different developmental stages (31). We found a significant overlap between genes activated during postnatal β-cell maturation and genes repressed in replicating β-cells. Similarly, genes repressed in maturing β-cells were significantly overrepresented among genes induced in replicating β-cells (Fig. 2). As expected, genes downregulated during β-cell maturation and induced in replicating β-cells were enriched for cell cycle–associated GO categories. Further, genes induced in β-cell maturation and repressed in replicating β-cells were enriched for β-cell functions (Fig. 2). Together, these results indicate that the transcriptome of replicating β-cells reflects a less mature state.

To assess whether reduced expression of genes encoding β-cell function occurs in human β-cells, we examined the published transcriptome of a recently described human β-cell line upon a shift from proliferation to quiescence, induced by Cre-mediated excision of T-antigen and human telomerase reverse transcriptase (hTERT) (32). Cell cycle exit of these cells was accompanied by activation of genes associated with β-cell function, which significantly overlapped with genes repressed in replicating mouse β-cells (Supplementary Fig. 2). As expected, cell cycle genes activated in replicating mouse β-cells were silenced in the human β-cell line after growth arrest. Thus, reduced expression of genes encoding β-cell function during β-cell replication is a feature conserved between mice and humans.

To address the underlying molecular mechanisms, we searched for information about the regulation of genes repressed in replicating β-cells. Enrichment analysis revealed that repressed genes were regulated by key β-cell transcription factors. Genes repressed in replicating β-cells significantly overlapped with genes with promoters that were bound by the NK6 homeobox 1 (Nkx6.1) or pancreatic and duodenal homeobox 1 (Pdx1) proteins and with genes that were repressed in Nkx6.1, Pdx1, V-Maf avian musculoaponeurotic fibrosarcoma oncogene homolog A (MafA), and HNF1 homeobox A (Hnf1a) mutant β-cells (33–36) (Fig. 3). Notwithstanding this observation, RNA-seq data and quantitative real-time PCR (qRT-PCR) analysis revealed that Nkx6.1, Mafa, and Pdx1 transcripts are not significantly reduced in replicating β-cells (not shown). Similarly, neuronal differentiation 1 (NeuroD1), an additional transcription factor regulating β-cell function, is not significantly repressed in replicating β-cells. Moreover, FACS quantification of Pdx1 and Nkx6.1 proteins by coimmunostaining for these proteins and Ki67 revealed that their levels were identical in replicating and quiescent β-cells (Fig. 3). These results suggest that functional β-cell genes are reduced in replicating cells via post translational modifications of the responsible transcription factors or, alternatively, that other mechanisms act to reduce expression of functional genes.

Figure 3

Targets of Nkx6.1, Pdx1, and MafA are repressed in replicating β-cells. A: Genes repressed in replicating β-cells are significantly enriched for targets of the β-cell transcription factors. The overrepresentation significance was computed using a statistical test based on hypergeometric distribution. B: FACS histograms of islet cells immunostained for insulin, Ki67, and Nkx6.1 or Pdx1. Cells were gated on insulin+ β-cells. Controls represent cells immunostained with isotype control. wt, wild type.

Global Upregulation of Gene Expression in Replicating β-Cells

We noticed that most repressed genes in replicating β-cells were decreased approximately twofold (Supplementary Fig. 3). We reasoned that the transition from quiescence to proliferation might be accompanied by a global doubling in mRNA content per cell, as observed in studies of cell lines (37–39). The normalization method used in differential expression analysis software assumes that the compared cell populations produce the same amount of RNA per cell. Under this assumption, a situation where global mRNA content increases but a subset of genes maintains constant expression levels might be falsely interpreted as downregulation of the unchanged gene subset. One solution to this problem is to normalize the RNA-seq data to the number of cells, using spike-in standards reflecting cell number (40). However, our use of very small cell populations (∼1,000 GFP+ β-cells per preparation) precluded the quantification of cell number prior to RNA isolation and the use of RNA spike-in standards.

To test the hypothesis that genes showing apparent stable expression levels were in fact upregulated in replicating β-cells, whereas genes that appear downregulated were in fact unchanged in quiescent and replicating cells, we used smRNA-FISH (27). Using this approach, we quantified the number of individual RNA transcripts from selected genes in individual β-cells that were either quiescent or replicating. Dissociated islet cells from CcnB1-GFP mice were mounted and fixed and hybridized with probes to specific genes as well as the cell cycle marker Ki67, in addition to imaging of GFP fluorescence. Analysis of transcript number per cell showed that genes with expression levels that appeared unchanged in our RNA-seq data (beta actin [Actb], histidyl-tRNA synthetase [Hars], and paired box 6 [Pax6]) were in fact increased in GFP+ replicating cells. Strikingly, RNA transcripts from “repressed” genes (Glut2/Slc2a2 and glutamate-ammonia ligase [Glul]) remained constant in dividing β-cells (Fig. 4). When we normalized the number of Glut2 or Glul transcripts to the number of Actb transcripts per cell, we observed a statistically significant decline in Glut2 and Glul expression in replicating GFP+ β-cells (Fig. 4). The significant approximately twofold increase in Actb transcript levels in replicating β-cells was confirmed in more than five independent smRNA-FISH experiments. Based on this observation, we used the Actb gene in qRT-PCR analysis as a gene standard. We selected from the RNA-seq data a panel of unchanged genes and genes that showed decreased expression in GFP+ replicating β-cells. Normalizing to Actb, qRT-PCR analysis confirmed that unchanged genes (Hars; poly[A]-binding protein, cytoplasmic 1 [Pabpc1]; and solute carrier family 38, member 4 [Slc38a4]) were regulated similarly to Actb (i.e., were present in higher levels in replicating cells). In contrast, a set of genes encoding β-cell function (Glul; Ero1lb; solute carrier family 30 [zinc transporter], member 8 [Slc30a8]; glucokinase [Gck]; Cpe; Scg5; and Slc2a2) that were seen as downregulated in RNA-seq were reduced approximately twofold in replicating β-cells compared with Actb (Fig. 4), implying that they were resistant to the general increase in RNA. Together our data reveal differential transcription dynamics in replicating β-cells in vivo: transcript levels of most genes are approximately doubled, while a subset of genes, enriched for β-cell functions, is refractory to this global upregulation. The mechanism underlying this bimodal transcriptional regulation is still unclear. However, it implies that postmitotic daughter cells will carry reduced amounts of RNAs encoding for proteins involved in β-cell mature state and function. These transcripts must be replenished postdivision to maintain the differentiated status of the β-cells after they exit cell cycle.

Figure 4

Quantification of RNA molecules by smRNA-FISH in quiescent (GFP−) and replicating (GFP+) β-cells. A: Representative image of smRNA-FISH for the detection and measurement of single mRNA molecules of Actb, Glut2, and Ki67, as indicated. Actb level is elevated in GFP+ (arrow) compared with GFP− cells (arrowhead) (267 vs. 156 transcripts, respectively), whereas Glut2 levels remain unchanged (60 vs. 61, respectively). Ki67 (marker for replicating cell) is expressed only in GFP+ cells. Scale bar is 10 μm. Image is a maximal projection of 12 optical sections spaced 0.3 μm. B: Quantification based on 15–20 optical sections spaced 0.3 μm in the Z direction. Mean numbers (with SEM) of RNA molecules per optical Z-section per cell (≥20 GFP+ and GFP− cells) for the five measured genes. Each dot represents measurements from a single cell. C: The ratios of Glut2/Actb and Glul/Actb transcripts were computed per cell. Statistical significance was assessed using the nonparametric Mann-Whitney test (B and C) (*P < 0.05, ***P < 0.0001). D: qRT-PCR analysis of the indicated mRNAs. Relative mRNA levels (±SE) were calculated after normalization to Actb (n = 3) (*P < 0.05, **P < 0.005, paired t test). NS, not significant.

No Evidence for Aerobic Glycolysis in Replicating β-Cells

A hallmark of cell proliferation is aerobic glycolysis (the Warburg effect), thought to be important for allocating precursors to the biosynthetic needs of proliferating cells (41). Aerobic glycolysis occurs in tumors and in rapidly proliferating lymphocytes (42,43). Surprisingly, our transcriptome analysis did not reveal evidence for a glycolytic shift in replicating β-cells. Several glycolytic genes were in fact significantly repressed (Fig. 5). Relative downregulation of Gck, Glut2, and AldoA was confirmed by qRT-PCR (Fig. 4 and not shown). Moreover, lactate dehydrogenase A (Ldha), a key determinant of aerobic glycolysis, and the gene encoding lactate transporter (solute carrier family 16, member 1; [Slc16a1/Mct1]), which are both specifically repressed in β-cells (to allow for accurate coupling of glucose sensing to ATP generation), were not activated (Supplementary Table 2). Similarly the expression levels of genes encoding hexokinases (Hk1, Hk2, and Hk3) remained low in replicating β-cells. Interestingly, Ldha is expressed at high levels in the chronically proliferating human β-cell line and is repressed upon the shift of cells to quiescence (32). To further examine cellular metabolism in replicating β-cells, we stained islet cells with TMRE, a fluorescent probe for monitoring mitochondrial membrane potential. As expected, high glucose enhanced the TMRE signal (Fig. 5). Strikingly, replicating (GFP+) islet cells had a stronger TMRE signal than quiescent (GFP−) cells, indicative of enhanced mitochondrial activity and implying active oxidative phosphorylation (Fig. 5). These data suggest that the Warburg effect is not implemented in homeostatic proliferation of β-cells. Interestingly, a recent study showed that quiescent hematopoietic stem cells entering the cell division cycle also have increased mitochondrial activity (44). Potentially, aerobic glycolysis is a feature of frequent cell division cycles, as in cancer and immune cell activation, and not of quiescent cells bursting into rare cell division cycles. Further work will establish whether and how the metabolism of replicating β-cells is adapted to meet the biosynthetic needs of one cell turning into two.

Figure 5

Metabolic changes in replicating β-cells. A: Glycolytic genes repressed in GFP+ (replicating β-cells). Graph shows the mean expression (±SE) in GFP+ β-cells relative to quiescent GFP− β-cells, as derived from the RNA-seq data. All GFP+ values represent significant fold changes compared with GFP− values, as computed by DESeq (FDR <0.05). B–E: Mitochondrial membrane potential (∆Ψm) measured by TMRE staining in dissociated islet cells after incubation in low (3 mmol/L) and high (16 mmol/L) glucose. TMRE intensity is shown for GFP− (B) and GFP+ (C) islet cells incubated at both glucose concentrations; TMRE staining is shown in GFP− and GFP+ islet cells after incubation in low (3 mmol/L) and high (16 mmol/L) glucose (D and E, respectively).

Discussion

The Cell Cycle Machinery of β-Cells

Understanding the molecular biology underlying β-cell proliferation is a major challenge for regenerative biology. β-Cell lines are not suitable models, since by definition they have defective cell cycle control, and rare replicating β-cells in vivo have been inaccessible for systematic analysis. Using the CcnB1-GFP cell cycle reporter, we generated a comprehensive transcriptome of β-cell replication in vivo. Genes upregulated in replicating β-cells included members of the canonical cell cycle machinery, as well as other genes previously shown to be associated with cell proliferation (Supplementary Tables 3-1 and 4). In addition, multiple upregulated genes were previously associated with cancer but not with proliferation of normal cells (Table 2 and Supplementary Table 3-2). Our data suggest that these genes are novel components of the endogenous cell cycle machinery that likely affect tumorigenesis via impacting the cell cycle.

Furthermore, 140 genes activated at least twofold (beyond the global amplification of expression) were genes annotated with functions unrelated to the cell cycle, genes annotated with unknown function, or simply unannotated cDNA or predicted genes (Supplementary Table 3-3). We anticipate that those strongly induced in replicating β-cells (Table 1) are bona fide novel players in the cell cycle machinery, at least in β-cells. For example, katanin p60 subunit A-like 2 (Katnal2), which is associated with microtubule organization (GO classification), might be involved in β-cell mitosis. Mitochondrial fission regulator 2 (Mtfr2), with expression (both at the RNA and protein levels) highly correlated with that of well-characterized cell cycle genes (Cdc25c, Cdt1, and Kif18b), is another interesting novel effector of β-cell proliferation.

Ongoing studies with the CcnB1-GFP system in additional tissues are expected to define the global and cell type–specific in vivo proliferation signatures of normal and cancerous cells.

Global Upregulation of Gene Expression and Exclusion of Genes Encoding Differentiated Function

Replicating β-cells had reduced levels of genes involved in β-cell maturation and function. In principle, this could reflect the biology of a specialized subpopulation of replication-prone, less differentiated β-cells. However, previous studies have found no evidence for such a situation and instead concluded that β-cell mass homeostasis relies on infrequent proliferation of differentiated β-cells, which all have the same potential to divide (3,5,45). Thus, it is most likely that our data reflect the transcriptional changes that take place when a fully differentiated β-cell enters one cell-division cycle.

Using smRNA-FISH, we found that, in fact, “repressed” genes in replicating β-cells are not downregulated; rather, their mRNA levels remain constant in the face of a global doubling of mRNA levels. A global increase in mRNA levels in replicating cells has previously been described (37–39) and is thought to set the stage for the generation of two daughter cells with mRNA content identical to that of the mother cell. However, the exclusion from this amplification of genes encoding differentiated functions is novel and will require in-depth investigation. One mechanism that could account for the lower mRNA levels of some of the genes is late replication of their genomic loci during S-phase. This mechanism contrasts with the established positive correlation between early replication and transcriptional activity (46). In addition, Raj and colleagues recently uncovered a generic mechanism of transcriptional dosage compensation that renders mRNA levels insensitive to the precise timing of DNA replication during S-phase (39). An alternative mechanism could be that transcription factors activating the relatively “repressed’ genes are limiting, thus preventing the increased transcription despite increased DNA dosage. We note, however, that the cell cycle–related global RNA increase has been reported in cells transiting from G0 to G1 prior to DNA replication (37). An additional possibility is that differentiated genes have a distinct chromatin structure that underlies the lagging phenomenon. However, expressed housekeeping genes and expressed tissue-specific genes are not known to differ in chromatin structure. Another potential determinant of the process is Myc, a key transcription factor participating in cell replication and proposed to be a global amplifier of transcription in activated lymphocytes and embryonic stem cells (47,48). The seemingly repressed genes in replicating β-cells may be low-affinity Myc targets that are “left behind” during global amplification. Future experiments will attempt to distinguish between these and other models to understand how specific β-cell genes escape global transcriptional amplification during the cell cycle.

Regardless of the mechanism, our findings imply that a β-cell that has completed mitosis will have reduced levels of transcripts required for its function. The dynamics accounting for upregulation of these genes postmitosis are an interesting topic for further investigation. Finally, while replicating β-cells clearly maintain their overall identity and are not dedifferentiating, our results do show a fundamental distinction between the programs of differentiation and replication, as we have shown before for replicating hepatocytes (19).

The list of genes “repressed” in replicating β-cells is highly enriched for components of the physiological function of β-cells, namely, GSIS. We speculate that additional genes in this list, not previously associated with β-cell function, are a rich source of novel candidate regulators of β-cell function. Ongoing work in our laboratory is testing this idea systematically.

In conclusion, we describe a comprehensive catalog of the transcriptome of replicating β-cells in vivo, including multiple upregulated genes not associated previously with β-cell replication. This resource will be useful for future investigations of determinants of β-cell replication toward the identification of targets for boosting β-cell regeneration in diabetes.