Classifying Biophysical Subpopulations of Insulin Secretory Granules using Quantitative Whole Cell Structure Analysis

SXT provides multiparameter descriptions of ISGs

We first summarize and highlight the ISG parameters that can be extracted using SXT. With the full rotation SXT method, we were able to image the entire volume of a β-cell (Figure 1A). SXT tomograms have an isotropic resolution between 30-35 nm. The LAC-based contrast mechanism of SXT allows us to identify a range of ISG densities, which we interpret to be a range of ISG maturities (Figure 1B).³² Tomogram segmentations resulted in label fields of key organelles such as the ISGs, mitochondria, nucleus, and plasma membrane (PM) across whole cells (Figure 1C).

Individual ISGs were parameterized by their LAC values, size, and distances from other important organelles (Figure S1A). Each ISG has a distribution of LAC values (Figure 1D), which allows for the calculation of a LAC mean, standard deviation, skew and kurtosis (Figure S1B). Skew describes the asymmetry of the LAC distribution, indicating whether the granule contains proportionally more high- or low-density voxels, while kurtosis quantifies the peakedness of the distribution, reflecting how sharply concentrated or broadly dispersed the LAC values are within each ISG. To account for the ellipsoid shape of some granules, the geometric mean diameter was used to represent ISG diameters. We computed comprehensive distances between ISGs and cellular structures like the mitochondria and PM by using the Euclidean distance transformation (EDT) of the corresponding organelle masks (Figure 1E, F, S1C). Distance from the PM is a key parameter, since it provides insight into the mobilization of ISGs for secretion. To account for the relative shapes and sizes of the PM and nucleus found within cells, distance the PM was normalized to the inverse EDT of the nucleus (Figure S2). Due to pressure caused by the capillary plunge-freezing process, some of the cells have an indentation or roughly triangular shape (Figure 1A). A potential artifact of this process is that ISGs in regions that are “squeezed” may be displaced disproportionately compared to granules in other cellular regions, thereby affecting their relative location in the normalized PM EDT. These approaches collectively yield a set of parameters used to identify unique ISG subtypes: LAC mean, standard deviation, skew, kurtosis, diameter, and distances from the PM and mitochondria (Figure S3).

Clustering framework to define ISG subpopulations

After establishing a list of parameters for each ISG, we set out to create a comprehensive unsupervised clustering methodology to classify ISG subpopulations. We reasoned that groupings of ISGs should remain consistent even if the underlying data is perturbed in a non-essential manner such as bootstrapping.³³ To implement this reasoning, cluster stability analysis was used. This analysis validates that over many bootstrapped datasets, the cluster assignments for individual ISGs are relatively stable as measured by the Jaccard coefficient.³⁴ Additionally, to create generalizable ISG subpopulations, we used multiple cells to average over individual cell variation and responses to experimental conditions. In total, we analyzed ~36,000 ISGs from 44 INS-1E cells in seven experimental conditions from the Deshmukh et. al. dataset (Table S1). We aggregated ISGs from all the cells across the seven conditions and performed clustering on this pooled dataset to create our final subpopulation assignments. In this approach, subpopulations are composed of ISGs with similar characteristics across multiple INS-1E cells.³⁵

After applying stability analysis and k-medians clustering on the aggregated dataset, we found 5 ISG subpopulations (Figure S4A). ISGs were visualized and labeled by their assigned subpopulation within multidimensional uniform manifold approximation and projection (UMAP) space (Figure 2A). When ISGs are colored by their individual parameter values such as LAC mean, UMAP illustrations enable association between specific ISG features and subpopulation identity (Figure S4B). To verify the UMAP subpopulation structure, we also performed principal component analysis (PCA). When ISGs are colored by subpopulation identity in the PCA projection, distinct clusters are still visible (Figure S5A). Granule-intrinsic characteristics such as LAC mean, LAC standard deviation, and diameter account for 30% of the explained variance, while ISG location characteristics contribute an additional 18% (Figure S5B-C, Table S2).

Subpopulations were labeled in order of increasing average LAC mean to reflect relative ISG maturity (Figure 2B, Video S1). Additionally, one potential limitation of grouping ISGs across cells is the risk of misrepresenting individual cells. To validate the applicability of this approach, we found that these subpopulations broadly capture the characteristics of each INS-1E cell in our dataset (Figure S6).

We interpret Subpopulation 1 as an immature pool of ISGs primarily due to their reduced LAC mean (Figure 2B). Although ISGs with minimal biochemical density may go undetected using SXT, integrated structural models of granules at earlier maturation stages suggest they still possess appreciable LAC content.³² Subpopulation 2 ISGs are close to the plasma membrane and compared to granules in Subpopulation 1 have relatively higher LAC and diameter. A similar pattern was observed in a previous cryo-electron tomography (cryo-ET) study, where the authors found a greater proportion of larger, mature ISGs near the cell periphery.³⁶ Therefore, these granules are likely mature ISGs with a PM localization comparable to the readily releasable pool. Subpopulation 3 is far from the PM, so this subpopulation may represent the internal reserve pool of ISGs.³⁷ While Subpopulation 4 ISGs have smaller diameters than Subpopulation 3 ISGs, they exhibit a higher LAC mean and reduced LAC standard deviation phenotype. This suggests that these are condensed granules with a large core size relative to total granule size. The last grouping, Subpopulation 5, is notable since it is far from the plasma membrane yet contains ISGs with very high LAC mean, larger diameter, and elevated LAC standard deviation (Figure 2B,C). Increasing the display range of maximum LAC illustrates the high LAC standard deviation of Subpopulation 5 ISGs, indicating that these granules occupy a unique state in the ISG maturation process (Figure S7). As a result, we consider this subpopulation to be a secondary, smaller group (9% of total ISGs) of mature and highly dense ISGs held in reserve.³ While it is possible that lysosomes could be mislabeled as large ISGs, prior correlative cryo-fluorescence analysis indicates that potentially mislabeled lysosomes would make up a negligible proportion of the overall dataset.²⁴

β-Cell stimulation conditions remodel ISG subpopulation distributions

To understand biologically meaningful changes in ISG subpopulations, we analyzed ISG subpopulations across six stimulation conditions associated with different cellular signaling pathways along with no stimulation as a control (Figure S8). Grouping by experimental condition reveals that different insulin secretagogues shift the relative distribution of ISG subpopulations (Figure 3A). This finding is visually apparent when ISGs in representative cells from each condition are colored by subpopulation identity (Figure 3B, Video S2).

In the unstimulated condition, the majority of ISGs are immature granules contained in Subpopulation 1. To stimulate INS-1E cells, the glucose condition was used because glucose treatment is the main mechanism for insulin secretion via glucose-stimulated insulin secretion (GSIS). While this condition was still majority Subpopulation 1 and not statistically significant compared to the unstimulated condition (Figure 3A), there was also a greater proportion of condensed and mature ISGs. Similar ISG subpopulation shifts were also observed using methods that separate organelles based on their electrokinetic properties.¹⁸ Glimepiride (GLM), in contrast, is an antidiabetic sulfonylurea that increases insulin release through K_ATP mediated membrane depolarization.³⁸ Due to the mechanism of action of sulfonylureas, we anticipated and observed their subpopulation composition was not statistically significant compared to the unstimulated condition. Insulin secretory stimuli that act on membrane depolarization lead to the secretion of a wide range of ISGs based on their pH profile compared to glucose stimulation.⁴ This could explain the slight differences in subpopulation distribution between these two conditions.

Using glucokinase activator (GKA50), which increases the metabolic activity of glucokinase resulting in increased ATP production,³⁹ we observed a statistically significant increase in the condensed ISG subtype Subpopulation 4 and the highly mature subtype Subpopulation 5. In this condition, the main effect was an increase in dense core formation in both large and small granules. TAK-875 is a G-protein coupled receptor 40 (GPR40) agonist acting on free fatty acid receptor signaling pathways like inositol triphosphate (IP3) and diacylglycerol (DAG) signaling.⁴⁰ TAK stimulation resulted in a statistically significant shift to the mature subtypes, Subpopulation 4 and 5. Secretion of Subpopulation 5, a highly mature ISG subtype, would help explain the observed effect of augmented insulin release during the second phase of prolonged insulin secretion upon TAK administration.^41,42

Incretins are a class of gut hormones responsible for decreasing blood glucose levels by increasing insulin secretion. Here, we use Gastric Inhibitory Polypeptide (GIP) and Exendin-4 (Ex-4, a glucagon-like peptide-1 receptor (GLP-1R) agonist), which both increase insulin secretion via cyclic-adenosine monophosphate (cAMP) signaling.⁴³ Despite their similar signaling pathways, the subpopulation profiles of GIP and Ex-4 are different. The main distinction is between the relative ratios of Subpopulations 1 and 4. Subpopulation 1 is nearly absent and Subpopulation 4 is notably enriched in the Ex-4 condition, while Subpopulation 1 is still a major component in the GIP treated cells. The intracellular insulin to proinsulin ratio is higher in Ex-4 treated INS-1E cells relative to GIP treated cells.²⁴ Therefore, Ex-4 could significantly enhance the ISG maturation conversion process in a distinct manner from other drug stimuli.

Heterogeneity of Docked ISGs

The readily releasable pool of ISGs is often interpreted to be the subset of ISGs primed and docked for release at the plasma membrane (Figure 4A,B, Video S3).^2,44,45 In this context, “docked” granules are attached to the plasma membrane via the SNARE complex.⁴⁶ We sought to understand if there was meaningful heterogeneity in this docked pool of ISGs beyond its exocytotic definition by employing a clustering strategy.⁴⁷ To combine the mixture of continuous information such as LAC mean and categorical information such as docked or undocked, we used the Gower's distance metric (Figure S9A).^48,49Here, we used a threshold of ≤ 50 nm to define docked granules. At the resolution of these SXT tomograms (30–35 nm), a distance threshold of 50 nm represents a conservative estimate for capturing a range of ISG states associated with the plasma membrane. An important caveat in this definition is that at this threshold, we cannot rule out if the ISGs in this region are ready for release, involved in kiss and run secretion, or tethered to the membrane.²

After performing k-medians clustering on the ~36,000 ISG dataset using Gower’s distance and focusing on inherent features of ISGs, we found three classes of docked ISGs (Figure 4C, Figure S9B).⁵⁰ While our earlier analysis sought to establish an appropriate level of ISG heterogeneity, here we instead focus on optimizing the relative weighting of the binary docked feature in the clustering framework. We find a straightforward interpretation of the docked classes, with an increase in ISG LAC mean, LAC standard deviation, and granule diameter from Classes 1 to 3. Over 46% of all docked ISGs were in Docked Class 1, reflecting the fact that most docked ISGs originate from conditions where the INS-1E cells are treated with stimuli that do not increase maturity. Similar to Subpopulation 5 from the previous analysis (Figure 2B), Docked Class 3 is distinguished by elevated LAC mean, diameter, and LAC standard deviation.

When docked ISG pools are categorized according to experimental condition, distinct profiles emerge (Figure 4D). The composition of the unstimulated and glucose condition profiles were similar to each other, with the majority of both being composed of Docked Class 1. This aligns with our expectation for the unstimulated condition, where many ISGs would be primed for release at the PM. Due to the membrane depolarization mechanism of GLM, very few docked ISGs were found. The GKA treated cells have few docked ISGs, suggesting that this stimuli mechanism may alter the dynamics of secretion and replenishment of newly docked ISGs at a 30 minute time point relative to glucose stimulation alone. The GIP condition was primarily composed of Docked Class 1, potentially reflecting an increase in new ISGs that are trafficked to the cellular periphery. In contrast, the TAK and Ex-4 stimulation conditions result in an increase in Docked Classes 2 and 3. While ISGs from TAK stimulated cells are relatively large, increasing the amount of potential surface area that can be in contact with the PM, ISGs from the Exendin-4 condition with generally smaller diameters still display increased PM contacts. Therefore, these findings are not solely due to increased ISG size and point to more dense, mature ISGs ready for release in the second phase of insulin secretion.

ISG Subpopulations in Primary β-Cells

To evaluate the ISG subpopulation clustering process in a more biologically relevant model, ISGs from mouse primary β-cells were used. Notably, INS-1E cells contain less ISGs and insulin content on average compared to normal β cells.²¹ When we repeated the analysis on six unstimulated and two Ex-4 stimulated β-cells, totaling to ~18,000 ISGs, we found 5 clusters using clustering stability analysis (Figure S10A). The distribution and values of unnormalized parameters was similar between the INS-1E and primary cells to allow direct comparison of subpopulations (Figure S11). ISG Subpopulations 2, 3, and 5 in INS-1E cells (Figure 2B) appear to correspond to Primary Subpopulations 2, 3, and 5 in β-cells based on similarities between the z-scores of the five primary ISG parameters (Figure S10A). In contrast, Primary Subpopulation 1 ISGs have larger diameters and higher LAC standard deviation compared toINS-1E Subpopulation 1 but have comparatively similar LAC means. These variations in biophysical ISG classification may reflect the differences in GSIS response between primary β-cells and insulinoma INS-1E cells.⁵¹

When primary ISG subpopulations are segregated by condition (Figure S10B), we observe subpopulation shifts that help validate results from the INS-1E cells (Figure 3A). While not statistically significant, there is a reduction in the proportion of lower LAC mean ISG subpopulations following Ex-4 stimulation. Additionally, there is an increase in Primary Subpopulation 3, which are smaller, condensed ISGs. These observations can be visually demonstrated with cellular renderings of an unstimulated and Ex-4 stimulated primary cell (Figure S10C). Overall, this data displays similar trends to the INS-1E cells but suggests there may be slight variations in biochemical content between the two cell types.

Volumetric EM facilitates a complementary perspective of ISG subpopulations.

To further test the generalizability of this analysis pipeline, we also analyzed a vEM dataset. While this method lacks the quantitative biochemical density information of SXT, more subcellular structures can be visualized in vEM due to the resolution and membrane contrast principle of this method (Figure 5A). An important metric captured by vEM is the distance of ISGs from the Golgi apparatus, which plays a central role in ISG biogenesis and subsequent maturation.⁵² Additionally, completely immature ISGs not easily detected by SXT may be visualized with EM. In Müller et al., the authors used 3D FIB-SEM to map seven mouse primary β cells (3 low glucose, 4 high glucose for 1 hour) and described a variety of subcellular organelle rearrangements.³⁰ Briefly, primary mice islets were frozen with a high-pressure freezer and freeze substituted using a cocktail containing osmium tetroxide and uranylacetate to visualize cellular membranes. After Durcapan embedding and coating, the EM-prepared islet samples were subjected to repeated rounds of SEM imaging and Ga+ surface ablations of 4nm. Voxel dimensions were rescaled, resulting in an isotropic resolution of 4nm. Individual islet samples were imaged over the course of 2 weeks, allowing the imaging of large regions spanning tens of microns and the collection of whole volume β-cell datasets. We set out to apply clustering methods on this vEM dataset to identify subpopulations of ISGs that provide a complimentary perspective to the subpopulations found through SXT. In EM, ISGs can be characterized by diameter as well as by their distance from the PM, mitochondria, Golgi, and microtubules. The data parameterization and clustering stability process for continuous parameters was repeated, and 5 stable EM subpopulations were found (Figure 5B).

Comparison of EM Subpopulations with SXT Subpopulations from primary β-cells (Figure S10A, Table S3) adds context that neither method could achieve alone. To perform this analysis, we indirectly compare subpopulations based on analogous parameters from different cells. A summary of the ISG subpopulations from different imaging modalities can be found in Table S3. EM Subpopulation 4 are medium sized ISGs near the Golgi apparatus and microtubules, likely representing newly synthesized ISGs. This grouping corresponds to Primary SXT Subpopulation 1, defined by low LAC, linking newly formed ISGs to immature dense core granules. Primary SXT Subpopulation 2 and EM Subpopulation 1 can also be matched together, as both are near the PM and have average-sized diameters. These characteristics suggest that they are readily-releasable pool ISGs. EM Subpopulation 2 and Primary SXT Subpopulation 5 both consist of large ISGs. Their extended distance from the Golgi and elevated LAC values, respectively, suggest a population of highly mature, dense ISGs.

Additionally, statistically significant shifts in EM Subpopulations are apparent between two glucose-based conditions (Figure 5C). Subpopulations 4 and 5 increase in the high glucose condition, which indicates newly synthesized ISGs being created at the Golgi and trafficked along microtubules. Interestingly, the high glucose condition has a lower proportion of Subpopulation 2 compared to low glucose. This reduction may reflect the depletion of mature ISGs through insulin secretion under prolonged stimulation. Although this trend differs from our SXT findings, it may be explained by the difference in treatment duration, as these cells were exposed to glucose for 1 hour compared to 30 minutes in the SXT experiments.

Limitations of the study

When using SXT to image INS-1E and primary β-cells, inherent constraints in spatial resolution and contrast likely limits the detection of ISGs, especially immature populations.³² For instance, β-cells imaged using vEM tend to have many more ISGs per cell compared to β-cells visualized with SXT.³⁰ Therefore, when unsupervised clustering is performed on ISGs from SXT-imaged cells, the subpopulations may be biased towards mature ISG states. While ISGs obtained from vEM data may be more representative of the immature ISG population, the distinction between immature and mature granules is less clear without the LAC parameter. To improve the description of ISG maturity in EM-based datasets for subpopulation analysis, a categorical classifier based on ISG core shapes could be used.

While whole cell, mesoscale datasets provide a plethora of information per cell, the overall number of cells is limited in this study. In total 44 INS-1E (SXT), 8 primary β-cells (SXT), and 7 primary β-cells (vEM) for a total of 59 whole cells were analyzed. While sufficient to describe several ISG states, imaging and processing more cells per condition could validate our findings and potentially enable us to discover new ISG subpopulations with novel biological roles. Advances in vEM image acquisition speed and automated segmentation using machine learning could increase throughput in large-volume EM, while increased accessibility of SXT technology could increase the sample size of cells to study.^60,61

METHOD DETAILS

QUANTIFICATION AND STATISTICAL ANALYSIS

Statistical analysis was conducted using R version 4.4.1. To perform statistical testing on ISG subpopulation shifts, permutational analysis of variance (PERMANOVA) testing using the packages vegan and pairwiseAdonis were used. Homogeneity of Multivariate Dispersion checks with TukeyHSD were first applied to understand differences in within-group dispersion. Then, pairwise PERMANOVA was applied to determine statistically significant ISG subpopulation composition differences between the conditions. Statistical significance was set at *p < 0.5, **p < 0.01, *** p < 0.001 for all analyses.