Microscopic image-based covariation network analysis for actin scaffold-modified insulin signaling.

To infer a "live" protein network in single cells, we developed a novel Protein Localization and Modification-based Covariation Network (PLOM-CON) analysis method using a large set of quantitative data on the abundance (quantity), post-translational modification state (quality), and localization/morphological information of target proteins from microscope immunostained images. The generated network exhibited synchronized time-dependent behaviors of the target proteins to visualize how a live protein network develops or changes in cells under specific experimental conditions. As a proof of concept for PLOM-CON analysis, we applied this method to elucidate the role of actin scaffolds, in which actin fibers and signaling molecules accumulate and form membrane-associated protein condensates, in insulin signaling in rat hepatoma cells. We found that the actin scaffold in cells may function as a platform for glycogenesis and protein synthesis upon insulin stimulation.

Introduction

Extracellular signaling is a critical process for cell-cell communication in multicellular organisms. Following initial binding and receptor activation of the signaling molecule, specific signal transduction proteins downstream of the receptors are sequentially activated/inactivated in a cascade response, resulting in real-time changes in their abundance, state (e.g., post-translational modifications such as phosphorylation and SUMOylation), and localization in signal-recipient cells. In the case of insulin signal transduction, drastic morphological changes in both the plasma membrane and cytoskeleton (primarily actin filaments) are induced simultaneously, and function in the modulation of insulin signaling. For example, PI3K-mediated actin remodeling upon insulin stimulus is reportedly coupled with the glycolysis by releasing the actin-bound aldolase (Hu et al., 2016). Furthermore, there is a dramatic change in protein localization in insulin-stimulated cells. Therefore, both protein localization and information on post-translational modification state are as important as protein expression level. However, such valuable post-translational modification and localization (spatial) information remains to be actively determined and used toward elucidating the protein networks of signal transduction.

Immunofluorescence (IF) data of individual cells is likely to be one such dataset that contains such information. Using an IF method with validated antibodies, the Human Protein Atlas project revealed that many human proteins localize to multiple compartments rather than a single compartment in individual cells, indicating that signal transduction should be elucidated on the basis of information on the spatiotemporal dynamics of proteins in a single cell (Thul et al., 2017). Furthermore, novel image-based technologies, such as imaging mass cytometry, cyclic IF, and multiplex profiling of subcellular protein distribution, have revealed comprehensive cell phenotypes and subcellular spatial markers using more than 40 parameters derived from immunostained proteins obtained from single cells (Giesen et al., 2014; Lin et al., 2015, 2018; Bodenmiller, 2016; Goltsev et al., 2018; Lundberg and Borner, 2019; Schürch et al., 2020). These elegant single-cell analytical methods are suitable for clustering cell populations at a specific moment in time by visualizing the high dimensional data obtained from the multiple staining of different proteins at the single-cell level. Admittedly, especially in signal transduction, temporal differences in the abundance (quantity), modification state (quality), and localization of proteins upon stimulus are crucial toward determining specific responses of individual heterologous cell types. However, these temporal aspects have not garnered the same attention as image-based cell analysis.

Recent advances in microscopic analysis and improvements in computing power have allowed for the feasible extraction of large sets of quantitative data concerning the quantity, quality, and localization information of target proteins from microscope immunostained images of individual cells. Here, we derived a novel protein network linked to protein localization and modification information, termed Protein Localization and Modification-based Covariation Network (PLOM-CON). This network, also referred to as a covariation network, is based on the concept that one protein group, whose members are involved in the same signal transduction pathways, exhibits synchronized time-dependent dynamic behaviors at an optimal cellular localization, where the members function at their maximum.

In this study, we used PLOM-CON analysis and graph clustering to elucidate the biological function of the undefined actin scaffolds that transiently form in insulin-stimulated rat hepatoma H4IIEC3 cells. Our results show a novel function of the actin scaffold, acting as an intracellular platform for the initiation of Akt-dependent glycogenesis upon insulin stimulation through recruiting and localizing phosphorylated GSK3β (Ser9) and GYG2, both of which are involved in the initiation reaction of glycogen synthesis, into actin scaffolds. Consequently, our findings demonstrate that PLOM-CON analysis is a suitable method for the functional analysis of the transient localization-dependent protein functions in signal transduction in cells.

Results

Actin scaffold formation in hepatoma-derived H4IIEC3 cells upon insulin stimulation

We have previously shown frequent, prominent accumulation of round-shaped structures of actin fibers in hepatoma-derived H4IIEC3 cells upon insulin stimulation (Kano and Murata, 2019). Here, live imaging using H4IIEC3 cells transiently expressing LifeAct-EGFP (Riedl et al., 2008) revealed, upon insulin stimulation, the rearrangement of actin filaments and subsequent actin waves just beneath the plasma membrane (Figure 1A, Video S1). Using confocal microscopy, we found transient convergence of actin scaffolds in the cytoplasm (Figure 1A: arrows, Video S1). Actin scaffolding was observed in almost all cells at least once during image acquisition (∼60 min), although scaffolding formed most frequently 5-15 min after insulin treatment (Figure 1A, lower panel; Video S1). Further, IF revealed that a variety of proteins, which are reportedly activated by insulin stimulation, colocalized with actin scaffolds (Figure 1B), including active (phosphorylated) forms of the insulin receptor (pTyr999 and pTyr1185), Akt, phosphorylated Akt (pSer473), phosphorylated GSK3β (pSer9), mammalian target of rapamycin (mTOR), phosphorylated mTOR (pSer2448), phosphorylated mitogen-activated protein kinase kinase (MAPKK, also known as MEK) (pSer298), and phosphorylated extracellular-signal-regulated kinase (ERK) (pThy202 and Tyr204) (Figure 1B). Using high-resolution microscopy, we found actin mesh-like structures of the actin scaffold with Akt (pSer473) aligned on the actin fibers (Figure S1A). In a previous study (Kano and Murata, 2019), accumulation to the actin scaffold of actin-modulating proteins, such as sorting nexin 9 (SNX9) and actin-related protein 2 and 3 (Arp2/3), and endosomal proteins, such as early endosome antigen 1 (EEA1) and adaptor protein phosphotyrosine interacting with PH domain and leucine zipper 1 (APPL1), were also observed. These findings suggest that the actin scaffold is a transient membrane-associated structure, likely derived from multiple early endosomes surrounded by an actin mesh, created to accommodate the transient accumulation of a variety of signaling proteins to modulate insulin signal transduction. In addition, we found that exogenously added FITC-dextran was incorporated inside some actin scaffolds (Figure S1B). These results indicate that the actin scaffold may be comprised of concave (Ω-shaped) and invaginated membrane structures, likely several fused early endosomes, surrounded by an actin mesh (Figures S1A and S1B). Thus, we next used PLOM-CON analyses to elucidate the biological function of the actin scaffolds, which are difficult to purify and typically unsuitable for biochemical analyses.

Figure 1

Transient formation of actin domain structures upon insulin treatment in H4IIEC3 cells

(A) After transfection of the LifeAct-EGFP plasmid and serum starvation, live-cell time-lapse imaging of H4IIEC3 upon insulin treatment was performed. Images were acquired every 3 min after 1 μM insulin stimulation. Arrows indicate the actin structures. In the lower column, the period when the structures were observed in each cell are shown as arrows. Bar = 5 μm. (B) Serum-starved H4IIEC3 cells were stimulated with 1 μM insulin for 15 min, fixed, permeabilized, and stained with antibodies against Akt, phosphorylated Akt at T308 (Akt pT308), phosphorylated Akt at S473 (Akt pS473), phosphorylated GSK3β at S9 (GSK3β pS9), mTOR, phosphorylated mTOR at S2448 (mTOR pS2448), phosphorylated Erk at T202 and Y204 (Erk pT202Y204), EEA1, APPL1, phosphorylated insulin receptor at Y999 (IR pY999), or phosphorylated insulin receptor at Y1185 (IR pY1185) (green) and rhodamine-conjugated phalloidin (Actin, orange). Bar = 5 μm.

Video S1. Dynamics of LifeAct-EGFP in insulin-treated H4IIEC3 cells, related to Figure 1

LifeAct-EGFP-expressing H4IIEC3 cells were serum starved overnight, and stimulated with 1μM insulin. Time-lapse images were taken every 1 min using an A1+ confocal laser microscope system (Nikon). The movie is played back at 8 frames per second.

Experimental workflow of PLOM-CON analysis

The experimental workflow of PLOM-CON analysis is shown in Figure 2, and consists of several steps starting with sample preparation and antibody selection, as well as image acquisition and analysis, covariation network inference, and cluster generation and analysis. It is based on the concept that one protein group whose members are activated in the same signal transduction pathway exhibits synchronized time-dependent behaviors at an optimal cellular localization. PLOM-CON analysis creates a covariation network in which the nodes of proteins are connected when the multimodal spatiotemporal correlated behavior of proteins is detected.

Figure 2

The scheme of Protein Localization and Modification-based Covariation Network (PLOM-CON) analysis

Immunostained cell samples were prepared in 96-well plates. Each sample was fixed at a different time in the horizontal axis direction and stained by different antibodies in the vertical axis direction. Then, fluorescent images of the cells were acquired using a confocal, laser-scanning microscope equipped with autofocus equipment. Each element in the matrix corresponds to a state observed cell on the left top well (and ). We defined as the quantified value of the -th feature of the -th protein of the -th compartment in the -th cell in the matrix using NIS-Elements ver 4.4 (Nikon). For example, a feature value obtained from imaging includes nuclear brightness, nuclear volume, cytoplasmic brightness, and so on. We computed the medians of each column, then had the variables with the indices of features, compartments, and proteins. We obtained the large matrix by merging each small matrix, after removing the variables with fewer signals. The large matrix included the whole information of the features, compartments, and proteins for each time. The graphical lasso was applied to the matrix, and then the sparse precision matrix was derived. Here, if an element of the matrix is zero, the corresponding variables are conditionally independent. After computing the partial correlation matrix from the precision matrix, we took the maximum value of each block matrix (e.g., P1C1) corresponding to the unit of a compartment in a protein. To obtain a network consisting of proteins, we only focused on the elements that represented the relationships between different proteins. The matrix was then visualized as a graph that we called the ‘covariation network’. We then performed the graph clustering by using the OCG method to extract the clusters.

Sample preparation for PLOM-CON analysis: Immunostained samples were prepared as follows: H4IIEC3 cells were cultured semiconfluently in 96-well plates, incubated at 37°C for 1, 5, 10, 15, 20, 30, 40, 50, or 60 min after 1 μM insulin stimulation, and then fixed. Samples were then processed using indirect IF methods. The cells were stained by the antibodies against protein of interest, anti-GAPDH antibody as a cytoplasm marker, Hoechst 33,342 as a nuclear marker, and rhodamine-phalloidin as an actin marker. The use of antibodies against post-translational protein modifications, such as phosphorylation, ubiquitination, and SUMOylation, enabled us to obtain multiple quantitative data points concerning relevant qualitative changes in proteins of interest and their localization over time. For this purpose, we chose a series of antibodies related to insulin signal transduction and the regulation of actin dynamics [refer to Section (3) for details]. Furthermore, we performed pathway enrichment analysis using the gene expression datasets generated from H4IIEC3 cells treated with DMSO or CK666, which disrupt actin scaffolds (refer to the next section), in the presence or absence of insulin treatment. According to the functional classifications in Kyoto Encyclopedia of Genes and Genomes, the enriched pathway that is mainly involved is “Insulin resistance (rno04931),” which is related to glycogen synthesis, and gluconeogenesis via activation of the PI3K-Akt pathway following insulin stimulation. This step is important because performing pathway enrichment analysis ensures the plausible selection of proteins stained in PLOM-CON analysis.

Image acquisition by confocal microscopy and image analysis: Fluorescent images of cells were acquired using a confocal, laser-scanning microscope equipped with an autofocus. Sectioning in the z axis direction was performed for almost complete cell coverage (see STAR Methods). We obtained up to ∼2,000 cells through automated image acquisition of multiple points in each well. Because intracellular localization of a target protein is one of the most important pieces of information in PLOM-CON analysis, we subdivided the cell into three compartments: the nucleus, cytoplasm, and actin scaffold, using the images of DAPI, GAPDH, and phalloidin staining (Figure 2). We measured the mean fluorescent intensity of each protein of interest for each cell area and for each of the three intracellular compartment areas (See “Quantification of images” section in Methods and Table S1). We also obtained the number and area of aggregates of several proteins (Table S1). To ensure the quality of segmentation, we carefully optimized the segmentation algorithm and further validated by examining the images for all experimental conditions (See “Automated confocal well plate imaging” section in STAR Methods). In addition, we examined the extent to which segmentation quality affected the resulting network by creating low-quality segmentation data (see “Validation of segmentation quality” section in STAR Methods and Tables S2 and S3 for details). As a result, we found that: 1) it is important to avoid the inclusion of the background, 2) smaller segmentation of cytosol than the actual area would not significantly influence our covariation network estimation, and 3) the nucleus' segmentation would not have a big influence on our method.

All of these feature quantities were acquired as single-cell data. The median for each well was then calculated as a representative value for each experimental condition because distribution of the data was typically skewed because of cellular heterogeneity. Based on time-series changes in each feature quantity of the target proteins, a heatmap of the median for each feature quantity displays substantial changes in quantity, quality, and localization of the target proteins upon insulin treatment (Figure 3A).

Figure 3

Covariation networks for insulin-stimulated H4IIEC3 cells

(A) A heatmap displaying the median for feature quantities in H4IIEC3 cells at different time points after insulin treatment. Red or blue indicate that the value was increased or decreased compared to 0 min, respectively. A, B, C, D, or E in the column of item represents synthesis mean intensity, mean intensity, sum volume, sum count, or count, respectively. (B and C) Covariation networks for insulin-stimulated H4IIEC3 cells at the ρ of 0.91. Proteins stained with antibodies are represented as nodes and the edge between the pair of nodes indicates that the feature quantities of the two proteins are correlated. The location at which the feature quantity was measured is indicated by the colored region (B) or the subnode at each end of the edge (C) (red: cytosol, orange: actin scaffold, blue: nucleus). (D) The network of pAkt (Ser473) with its neighbor nodes. (E) The nodes for Akt (pSer473), p70S6K, Akt, FoxO1, MEK1 (pSer298), GSK3β (pSer9), and EEA1 with the subnodes are magnified from the covariation network of Figure 3C.

Selection and validation of antibodies in the PLOM-CON analysis: Antibodies used in this study (Table S4) were carefully validated by IF method. The antibodies with which immunofluorescent images showed sufficiently high signals (S) compared with the background (N) (i.e., where the S/N ratio was 5 or greater) were used for the analysis. We also validated whether the appropriate staining patterns for the protein were observed. Figure S2 presents changes of the mean level or ratio of fluorescence intensity between the nucleus and cytoplasm of all phosphorylated proteins analyzed in this study upon insulin treatment. These findings indicate that changes in the quality (phosphorylation) and localization of target proteins can be captured at the single-cell level. Despite these criteria, we suppose that the complete elimination of the non-specific or false staining of proteins might be quite difficult. Dealing with the problems, we adopted the sparse estimation with the graphical lasso because it is expected to mathematically remove the effect of non-specific or false staining (for details, refer to the next section). Supported by mathematical principle, our plan was to use as many antibodies as possible, which were pre-screened using the aforementioned strategy, to obtain as many time-series of feature quantities as possible.

Covariation network inference: We distinguished signal values from noise for all feature quantities derived from cell fluorescent images. To validate these values, we used the mean fluorescent intensity of GAPDH as a noise judgment criterion, as described in the Methods. For all variables that were not eliminated by the abovementioned criteria, a matrix was created with feature quantities in the horizontal direction and time in the vertical directions. This matrix was used as the basis for inferring the relationship between the variables.

To infer a biological network, a commonly used method is correlation analysis such as Pearson's correlation. However, these traditional approaches can yield spurious correlations that are unable to distinguish between direct and indirect associations. To avoid this issue, Gaussian graphical models (GGMs) have been proposed and widely accepted for use in various domains (Friedman et al., 2008; Ma et al., 2007). In formulating the GGMs, assuming that the data samples are independently generated from a multivariate Gaussian distribution, the network is represented as the inverse of the covariance matrix or adjacency matrix. In the current study, we used a graphical lasso as one of the techniques based on sparse GGMs (SGGMs) (Idé, 2009; Montastier et al., 2015; Varoquaux et al., 2012) considering feature quantities as random variables. This method is expected to provide at least three benefits to our analysis: (1) To solve the singular covariance inversion problem, (2) to obtain meaningful selected edges in the network, and (3) to utilize as many antibodies as possible. The benefit of (3) is based on the following scientific rationale. If some antibodies exhibit non-specific staining, it is possible that time-series changes in the feature quantity of the images using these antibodies are caused by random events. In this study, we calculated the precision matrix via L1 regularization, thereby eliminating any spurious correlations derived from non-specific staining signals. Therefore, even if the feature quantities for non-specific staining have correlations with other feature quantities, it can be expected that such a correlation with the feature quantity does not appear in the precision matrix.

We obtained a partial correlation matrix by normalizing the precision matrix with diagonal components. As protein localization is essential information, we expressed the relationship between proteins as a network that included the localization information. We selected the maximum value of the partial correlation values for each compartment pairs between two proteins, and drew a graph comprising nodes corresponding to the proteins in the compartments, where the edges indicate a relationship among the proteins in the compartments. We termed the graph a covariation network. Figure 3B shows a covariation network with nodes grouped by subcellular localizations: cytoplasm, nucleus, and actin domain. The inter-protein relationships across the intracellular compartments were easily visible. This covariation network was also expressed as shown in Figure 3C, highlighting the intracellular compartments that were shown as subnodes at each end of the edge. The subnodes were color-coded according to the intracellular compartment. These network diagrams emphasized the relationships between the proteins considering the protein localization.

The graphical lasso includes a regularization parameter (ρ). This parameter controls the sparsity of the network, which means the network can have a small number of edges at higher values of ρ. Theoretically, a network with lower ρ does not necessarily contain all the edges with higher ρ. We therefore examined in our data whether an edge, once disappeared, may re-appear in a higher ρ in our covariation network by changing the value of ρ from 0.71 to 0.93 (DMSO) by 0.02. No edges re-appeared at higher ρ in our case.

To determine the ρ value as a data-driven method, we first tried the cross-validation (Figure S3 and Table S5). However, at the ρ value determined by cross-validation, the network was too dense for us to extract the biologically significant information from the network (Figure S3 and Table S5). Therefore, we adopted a criterion based on biological relevance to determine the ρ value.

The insulin signaling pathway comprises successive and yet complicatedlycrossed reactions between proteins, and insulin receptor (IR) is considered the most upstream protein. IR firstly receives insulin at the plasma membrane, while FOXO1 is one of the most downstream proteins regulating gluconeogenic genes' transcription in the nucleus. To reveal the proteins' correlation “within” the insulin signaling, we adopted a regularization parameter with a lower limit of magnitude that eliminates correlations between proteins in the most distant positions within the insulin signaling cascade. If we adopted the maximum regularization parameter in which the correlation between IR and FoxO1 remained, we would be overestimating the correlation inside the insulin signaling system. Conversely, if we choose the minimum value of ρ for which the edges between the most distant proteins (IR and FoxO1 in our case) disappear, the correlation between relatively distant proteins in the signaling system would be set to zero, and we can detect the correlation between closely related proteins. Therefore, by adopting this criterion in this study, we could obtain the most appropriate sparsity network for our aim, to estimate the signaling system's internal correlation. This criterion's validity is supported because no direct biological interaction between IR and FoxO1 has been confirmed.

A covariation network of insulin-treated cells is shown in Figures 3B and 3C (ρ = 0.91) and Figures S4–S6, demonstrating different ρ values. In the covariation network, protein pairs with edges may show similar behaviors upon insulin stimulation and are likely to belong to a similar group in the insulin signal transduction pathway. Proteins in a specific compartment correlated with others in various compartments, indicating that the proteins functioned synchronously across different compartments. An edge in the covariation network indicates the existence of a positive or negative correlation between two nodes. A heatmap of feature quantities can indicate whether a given correlation is positive or negative (Figure 3A). As shown in Figure 3D, Akt (pSer473) had numerous edges and correlated primarily with Akt, GSK3β (pSer9), p70S6K, FoxO1, and EEA1. Observing the subnodes of each protein, Akt (pSer473) included the subnode of all compartments, and p70S6K primarily appeared in subnodes of the cytoplasm and nucleus (Figure 3E). Akt, FoxO1, and MEK1 (pSer298) were in the subnode of the nucleus, whereas GSK3β (pSer9) and EEA1were primarily connected to other nodes via actin scaffold subnodes (Figure 3E). These findings strongly suggest the possibility that biological functions in which GSK3β (pSer9) and EEA1 play essential roles depends on the actin scaffold.

Graph clustering to select clusters from the covariation network: As the covariation network at low ρ appears too complicated to readily extract the biologically significant group of proteins, we next performed graph clustering to statistically extract network characteristic features, listing groups of proteins with similar physiological correlations. We used the overlapping cluster generator (OCG) method (Becker et al., 2012) which allows one node to belong to multiple clusters. After extracting a cluster by OCG, we then computed the mean number of edges within each cluster (mean degree of the cluster) and the mean degree of the whole covariation network for each ρ to extract the characteristic cluster. In general, we found that the mean degree of each cluster generated by OCG varied widely (Figure 4A). Therefore, we next extracted those clusters with a mean degree higher than that of the entire network (for each ρ) as characteristic clusters (Figures 4A and S7–S13).

Figure 4

Graph clustering of the covariation network for insulin-stimulated H4IIEC3 cells

(A) Clusters were detected using OCG algorithm from the covariation network. The mean degree, i.e. mean number of edges that the node has to other nodes, in each cluster were plotted (blue). The red line indicates the mean degree in the covariation network at each ρ. (B) Graph clustering of the covariation network for insulin-stimulated H4IIEC3 cells (ρ of 0.79). (C) A cluster with mean degree larger than that of the covariation network at the ρ of 0.79. The subnode at each end of the edge indicates where the feature quantity was measured (red: cytosol, orange: actin scaffold, blue: nucleus). (D) Graph clustering of the covariation network for insulin-stimulated H4IIEC3 cells (ρ of 0.93). (E) A cluster with mean degree larger than that of the covariation network at the ρ of 0.93. The subnode at each end of the edge indicates where the feature quantity was measured (red: cytosol, orange: actin scaffold, blue: nucleus).

At a ρ value of 0.79, we found that nodes were grouped into several clusters (Figure 4B), and a characteristically large cluster was extracted (Figure 4C). This cluster consisted of Akt, its activated form Akt (pSer473), FoxO1, FoxO1 (pSer256), GSK3β (pSer9), GYS2, MEK1, and EEA1. As shown in Figures S7–S13, as the ρ value increased, the number of proteins in the cluster gradually decreased. Interestingly, the clusters shown in Figure 4C could be divided into two smaller clusters with high ρ values. At a ρ value of 0.93, we found that the number of clusters decreased (Figure 4D), and two smaller clusters (Cluster #1 and Cluster #2) were obtained (Figure 4E). Akt, Akt (pSer473), and EEA1 were common members of both clusters. FoxO1 and FoxO1 (pSer256) were also included in Cluster #1, whereas GSK3β (pSer9) was in Cluster #2 (Figure 4E). GSK3β (pSer9) was linked through the subnode of the actin scaffold, but FoxO1 and FoxO1 (pSer256) were linked through the subnode of the nucleus (Figure 4E). These results suggest the involvement of actin scaffolding primarily in GSK3β (pSer9)-mediated glycogenesis.

Biological evaluation of localization- and actin scaffold-dependent glycogen synthesis

Possible actin scaffold-dependent biological processes were determined by focusing on those proteins connected to other proteins via an actin scaffold subnode, namely, GSK3β (pSer9) and GYS2 (Figure 4C). Further, as GSK3β (pSer9) and GYS2 have important roles in glycogenesis (Cross et al., 1995; Orho et al., 1998), we hypothesized that inhibition of actin scaffold formation affects glycogen synthesis. To test this hypothesis, we measured glycogen synthesis in insulin-stimulated cells, in which actin scaffold formation was inhibited. We used an inhibitor of Arp2/3 CK666, which inhibits the branching of actin fibers (Hetrick et al., 2013). CK666 effectively inhibited actin scaffold formation at 50 μM without inducing other substantial morphological changes in the H4IIEC3 cells (Figure S14). Further, insulin-stimulated glycogen synthesis was inhibited by CK666 treatment (Figure 5A), regardless of the expression level or phosphorylation state of GSK3β at Ser9 (Figure 5B). Thus, we reasoned that inhibition may be because of a localization-specific function of GSK3β (i.e., actin or the actin scaffold becomes a platform to enhance glycogen synthesis by recruiting these proteins). To examine whether glycogen synthesis-related proteins may also accumulate with actin scaffolds, we extracted the candidate proteins by expanding the covariation network through the addition of nodes for the interacting proteins (see STAR Methods for detailed information, Figure 5C). In the cluster including both GSK3β (pSer9) and GYS2, the candidate proteins that operate in glycogen synthesis and interact with GYS2 were GYG2, protein phosphatase 1 regulatory subunit 3C (PPP1R3C), PPP1R3B, PPP1CA, and starch binding domain 1 (STBD1). Among these proteins, we observed that GYG2 was distributed throughout the cytoplasm and accumulated with actin scaffolds upon insulin stimulation, but the others did not (Figure 5D). GYG2 is a self-glycosylating protein that initiates glycogen synthesis in liver cells (Mu et al., 1997). Therefore, the specific accumulation of GSK3β (pSer9) (Figure 1B) and GYG2 at actin scaffolds indicates that the domain may promote the accumulation of proteins related to initiation of glycogen synthesis and other upstream signaling proteins, leading to effective upregulation of insulin-mediated glycogen synthesis. In contrast, inhibiting gluconeogenic gene expression of phosphoenolpyruvate carboxykinase 1 (PCK1) and glucose-6-phosphatase catalytic subunit (G6PC) upon insulin stimulation occurred normally in both groups (Figure 5E).

Figure 5

Biological evaluation of the covariation network for insulin-stimulated H4IIEC3 cells

(A) Glycogen assay in DMSO- (blue) or CK666- (red) treated cells. Six independent experiments were performed. ∗P < 0.05. (B) DMSO- or CK666-pre-treated H4IIEC3 cells were incubated with insulin for the indicated time. Western blotting was performed using antibodies against GSK3β (pSer9) and GSK3β. The mean and the standard deviation of the ratio of GSK3β (pSer9) to GSK3β are shown in the graph. Three independent experiments were performed. (C) The expanded covariation network containing the first neighbor interacting proteins. In the cluster that included both GSK3β (pSer9) and GYS2 (cluster), the proteins that interact GYS2 and function in glycogen pathway were extracted (node). (D) H4IIEC3 cells were stimulated with insulin for 15 min, and was subjected to immunofluorescence. GYG2, Myc-tagged PP1R3C, PPP1R3B, PPP1CA, and Myc-tagged STBD1 were co-stained with rhodamine-phalloidin (Actin). Bar = 10 μm. (E) Real-time PCR analysis of G6PC in DMSO- or CK666-treated cells. Three independent experiments were performed. ∗P < 0.05. ∗∗P < 0.01.

Comparative analysis to elucidate the biological function of actin scaffolds

Next, we inferred the covariation network for CK666-treated cells for comparison. The network should reflect the impact of actin scaffold disruption on global cell responses. As expected, the heatmap of the median for feature quantities (Figure 6A) and the shape of the covariation network in CK666-treated cells (Figures 6B, S15, and S16) were different from those of DMSO-treated control samples (Figure 3C). Most notably, we found that Akt (pSer473) was connected to mTOR (pSer2448) and eukaryotic translation initiation factor 4E (eIF4E) in CK666-treated samples (Figure 6C). These combinations of proteins, including Akt (pSer473), mTOR (pSer2448), and eIF4E, were also extracted by graph clustering of the covariation network in CK666-treated cells (Figures 6D, 6E, S17, and S18). Compared to the covariation network in CK666-treated cells, we found that most of the edges of phosphorylated mTOR (pSer2448) were connected via subnodes of the actin scaffold in the covariation network in DMSO-treated cells (Figures S4 and S5, ρ = 0.79 and 0.85). We observed a connection between mTOR (pSer2448) at the actin scaffold with MEK1 (pSer298), MEK1, EEA1, mTOR, GSK3β (pSer9), and FoxO1/2 in DMSO-treated cells (Figure S5, ρ = 0.85). These edges were lost by CK666 treatment, and an Akt (pSer473)/mTOR(pSer2448)/eIF4E cluster appeared (Figure 6E). Furthermore, several proteins related to protein synthesis, such as p70S6K (pSer424) and S6RP (pSer235/236), were also connected with other nodes via the actin scaffold in DMSO-treated cells (Figures S4 and S5, ρ = 0.79 and 0.85). These results indicate that mTOR signaling and its downstream pathway to protein synthesis may be modified when actin scaffold is disrupted.

Figure 6

PLOM-CON analysis for insulin-stimulated H4IIEC3 cells in the presence of CK666

H4IIEC3 cells were pretreated with 50 μM CK666 or with DMSO (control) for 20 min and then further treated with 1 μM insulin at 37°C for 1, 5, 10, 15, 20, 30, 40, 50, and 60 min. The cells were fixed and subjected to PLOM-CON analysis. (A) A heatmap displaying the median for each feature quantity for CK666-treated H4IIEC3 cells. Red or blue indicates that the value was increased or decreased compared to that at 0 min, respectively. A, B, C, D, or E in the column of item represents synthesis mean intensity, mean intensity, sum volume, sum count, or count, respectively. (B) Serum-starved H4IIEC3 cells were pretreated with 50 μM CK666 for 15 min and stimulated with insulin, and was applied to PLOM-CON analysis. A covariation network at the ρ of 0.86 is shown. The location at which the feature quantity was measured is indicated by the subnode at each end of the edge (red: cytosol, blue: nucleus). (C) The network of Akt (pSer473) with its neighbor nodes from the CK666 covariation network in B. (D) Clusters were detected using OGC algorithm from the covariation network for CK666-treated cells. The mean degrees in each cluster were plotted (blue). The red line indicates the mean degree in the covariation network for CK666-treated cells at each ρ. (E) A cluster with mean degree larger than that of the covariation network for CK666-treated cells at the ρ of 0.88. The subnode at each end of the edge indicates the location at which the feature quantity was measured (red: cytosol, orange: actin scaffold, blue: nucleus).

Although insulin-stimulated new protein synthesis was undetectable in cells of both groups using commercially available probes, likely because of a lack of sensitivity, we succeeded in detecting less-phosphorylated mTOR at Ser 2448 (Figure 7B). Interestingly, we observed no significant differences in the kinetics of Akt phosphorylation between DMSO- and CK666-treated cells (Figure 7A), suggesting that actin scaffold disruption modulates signal transduction downstream of Akt in the insulin-signaling pathway, but not the phosphorylation state of Akt itself. As the phosphorylated form of mTOR at Ser2448 is reportedly included in the mTORC1 complex (Copp et al., 2009), we tested the phosphorylation state of its best-characterized downstream proteins, S6 kinase (S6K) and eIF4E-binding protein (4E-BP). Using western blotting, 4E-BP phosphorylation at T37 showed no significant change, but phosphorylation of S6K at Ser424 was slightly decreased by CK666 treatment (Figures 7C and 7D). Furthermore, we found substantial inhibition of phosphorylation of S6 ribosomal protein (S6RP) (Figure 7E), which is a direct substrate of S6K. These results suggest that CK666 treatment may negatively affect protein synthesis through repression of insulin-induced phosphorylation (activation) of mTOR and S6RP. These findings further indicate the possibility that the localization of mTOR to actin scaffolds affects not only its phosphorylation/activation state, but also the signaling efficiency of downstream protein synthesis.

Figure 7

Graph clustering and the biological evaluation of covariation network of insulin-stimulated, CK666-treated H4IIEC3 cells

(A-E) Serum-starved H4IIEC3 cells were pre-treated with DMSO (blue) or 50 μM CK666 (red) for 20 min, further with 1 μM insulin for 0, 5, 15, 30, 60 min, and electrophoresed protein extracts subjected to western blotting using antibodies against Akt or Akt(pSer473) (A, left), mTOR or mTOR(pSer2448) (B, left), p70S6K or p70S6K (pSer235/236) (C), 4E-BP or 4E-BP (pThr37/46) (D), and S6RP or S6RP (pSer235 and pSer236) (E, left). The mean and the standard error of the ratio of phosphorylated protein to total protein are shown in the graph (A, B, and E, right). Three independent experiments were performed. ∗P < 0.05. (F) Proposed model for the function of the actin domain in the insulin-signaling pathway.

Discussion

An aim of this study was to establish and demonstrate the applicability of a versatile analytical method, PLOM-CON analysis, to investigate microdomain-mediated signal transduction in cells. As a proof of concept, we first applied PLOM-CON analysis to identify the function(s) of actin scaffolds, which were transiently formed in H4IIEC3 cells upon insulin stimulation, and to elucidate the aberrant impact of their disruption or dysfunction on the insulin signal transduction network.

There are several advantages to studying actin scaffold function in cells using PLOM-CON analysis. First, information regarding post-translationally modified proteins may be included in the network. Information of modified proteins (especially, phosphorylated form) enables us to infer a novel protein network that can include regulatory intermediates for signaling as nodes. Thus, the use of antibodies provides such unique characteristics to PLOM-CON analysis, compared to the mass spectrometry-based spatial proteomics (Christoforou et al., 2016; Itzhak et al., 2017; Lundberg and Borner, 2019). In regards to this point, imaging mass cytometry using specially modified antibodies with isotopically rare earth metals or DNA barcodes has become increasingly common for multiplex profiling of phosphorylated intermediates of signal transduction for single-cell analysis (Bodenmiller et al., 2012; Bendall et al., 2011; Lun et al., 2017; Lundberg and Borner, 2019). Our covariation network information, which contains temporal synchronicity of the proteins in the same signaling pathway, provides corroborating information to the multiplex profile of intermediates in single-cell analysis. Second, information on protein (co)localization may also be included in the feature quantities and shown as specific subnodes in the covariation network (Figures 3 and 6). In addition to protein localization information, morphometric information (e.g., number, area, circularity, and elongation) and associated aggregate structures of proteins are also obtainable from immunostained cell images. Such image-based quantitative information allows us to infer organelle- and/or domain-specific covariation networks. Thus, as shown in Figure 3, the synchronous behavior of proteins between the actin scaffold, nucleus, and cytosol can be compared under various experimental conditions. Third, as PLOM-CON analysis is based on cell image data, this method is applicable to heterogeneous cell populations by classifying or extracting a subpopulation of cells having certain properties. For example, cells under different conditions, such as healthy/pathogenic or undifferentiated/differentiated states, or in co-culture systems that contain many types of cells, can be sorted based on their immunostained cell images or by a machine-learning-based image recognition system. Interestingly, the network is also able to show correlations between proteins not only within a single cell type, but also between different types of cells. Thus, by taking extensive use of these advantages, PLOM-CON analysis complements quantity, quality, spatiotemporal, and morphological information of proteins in actin scaffolds, which are exceedingly difficult to isolate and preserve using biochemical cell preparations for microarray or other comprehensive analyses such as mass cytometry (Bodenmiller et al., 2012; Bendall et al., 2011; Lun et al., 2017).

Although the use of IF methods confers several advantages, in general, it should be noted that appropriate care should be taken when quantifying and evaluating fluorescence data of IF experiments. Regarding the issue of quantitativeness of IF data, in this study we calculated correlation efficiencies between the time series of the value in each feature quantity and compared them with their temporal synchronism (Figure 2). This approach allows us to compare fluorescence intensities between different proteins without the need to consider the absolute quantity of proteins. For antibody validation, the antibodies used for quantitative analysis of IF images were often validated both biochemically and morphologically (Lundberg and Borner, 2019). However, antibody validation is neither an easy nor a simple task. For instance, it would be extremely difficult to simply compare florescence data obtained from individual immunostained cells to that of Western blot data, which are obtained biochemically as an average value of a large number of cells. In addition, as IF signals are usually dependent on antibody accessibility to the protein epitope, they may be affected by conformational changes of proteins caused by experimental processes such as fixation, interactions with other proteins, or oligomerization state. In other words, IF images may include a far greater and more valuable variety of molecular information than one may imagine; for example, information on the intracellular structural or environmental status of antigen proteins related to a signal transduction pathway. Therefore, in this proof of concept, we proactively took advantage of IF image data and used as many available antibodies as possible.

To do so, we used sparse estimation via L1 regularization to obtain the precision matrix (Figure 2). The sparse estimation method is designed to eliminate spurious correlations including the various inappropriate or confounding effects of non-specific antibody staining signals (see (3) Selection and validation of antibodies in the PLOM-CON analysis and (4) Covariation network inference in Results section). Based on the intended robustness of these underlying methods, we selected antibodies for IF that display a sufficient fluorescence signal in their respective immunostained image compared with background and exhibit appropriate subcellular staining for the protein of interest. As a result, respective covariation networks of insulin signal transduction and its perturbation by CK666 in hepatic H4IIEC3 cells were generated and several protein clusters of biological significance were extracted by graph clustering (Figures 4 and 6). Consequently, we found experimentally that the specific accumulation of GSK3β (pSer9) (Figure 1B) and GYG2 (Figure 5D) at actin scaffolds leads to effective upregulation of insulin-mediated glycogen synthesis. Based on these findings, our proof-of-concept study demonstrated feasible results, and therefore validated our antibody selection strategy.

We determined the minimum ρ based on biological relevance in this analysis. In some cases, it may be better to adopt the upper limit regularization parameter that correlates the most upstream and the most downstream. For example, when observing multiple signal transduction systems simultaneously, we are interested in the relationship between signaling pathways, such as crosstalk between signaling pathways, rather than focusing on a single signaling pathway and observing the details within that pathway. In such cases, it is reasonable to adopt the maximum ρ for detecting correlations between the most distant proteins in the observed signaling pathway so that correlations between signaling pathways can be observed. So, at this time, we believe that determining ρ values based on biological relevance is reasonable and that the most appropriate ρ value would depend on the study's purpose.

The PLOM-CON analysis system consists of two major processes: a data production process and a subsequent network inference process. In this study, conventional IF was used to acquire single-cell images in the former process, and graphical lasso or OGC was used in the latter. Importantly, imaging data obtained from eminent multistaining methods, namely, imaging mass cytometry (Bodenmiller et al., 2012; Bendall et al., 2011; Lun et al., 2017), cyclic IF (Lin et al., 2015, 2018), and modified DNA barcoding approaches (Goltsev et al., 2018; Schürch et al., 2020), are applicable for the data production process in our PLOM-CON analysis system, illustrating both the flexibility of the system and the broad scope of its potential use. Furthermore, the feature quantities obtained from the live imaging could be used in our PLOM-CON analysis, greatly improving the time resolution. MacHacek et al. (2009) performed live imaging of two GTPase biosensors coupled with a computational multiplexing approach, which revealed the precise spatiotemporal coordination of membrane protrusion and activation of actin dynamics-related three GTPases. So obtaining data from more time points with shorter time intervals would allow us to discuss the temporal order of protein activation/localization. This approach, however, includes several difficulties. For example, we have to establish the cell line in which the GFP sequence is inserted at the gene of interest to detect the change in its endogenous expression level. In addition, multiple cell lines that expressed GFP protein of interest would be crucial for obtaining the directed protein network in our study. Such experimental design is possible but probably laborious and time-consuming at this stage. Furthermore, detecting the intracellular behaviors of proteins having post-translational modification is one of the prominent features of PLOM-CON analysis. However, generating a GFP sensor for the modification (e.g., phosphorylation) of protein at the specific amino acid would be laborious. Another approach for improving time resolution is adapting a robot manipulating system for cell treatment and immunofluorescence. We performed all the experimental processes by hand (manually) in this study, but the robot manipulation system would reduce the manual-based processes and increase the time resolution of each sample preparation of immunofluorescence, at least at 30–60 s.

Regarding the network inference process, as the number of nodes is dependent on the number of antibodies used in each experiment, it is a fact that our covariation network lacks comprehensiveness (for details, refer to the next section). However, we posit that correlation analysis would reveal a synchronized link between the time-dependent behavior of proteins and that of other molecules such as metabolites and mRNAs, which are obtained by more comprehensive data acquisition. In addition, connecting the covariation network with a protein-protein interaction (PPI) database as we did in this study (Figure 5C) demonstrates an approach that would further add to the comprehensiveness offered by our PLOM-CON analysis system. Thus, the covariation network created by our PLOM-CON analysis has the remarkable capacity to unveil novel correlations not only between proteins in cells under specific conditions, but also between proteins and other molecules.

Furthermore, our new findings of the biological function of actin scaffolds in insulin-stimulated H4IIEC3 cells demonstrate a proof of concept for PLOM-CON analysis enabling the evaluation of localization-specific functions of proteins (Figure 7F). In fact, subnodes illustrated the localization-specific relationship between two connected proteins, and ultimately, the localization-specific function of different proteins of interest, including GSKβ (pSer9), GYG2, and mTOR (pSer2448). In this study, we first hypothesized from the subnodes of GSK3β (pSer9) and GYS2 that glycogen synthesis may be actin-scaffold-dependent. This hypothesis was verified by the following observation: (1) the inhibition of actin scaffold-formation attenuated insulin-mediated glycogen synthesis (Figure 5A), and (2) morphological observations revealed the accumulation of GSK3β (pSer9) and GYG2, which primes the polymerization of glucose for glycogen synthesis, to actin scaffolds without changing the amount of GSK3β (pSer9) (Figures 1B, 5B, and 5D). Indeed, at the onset of insulin-dependent glycogen synthesis, GYG2 binding to actin was reported in C2C12, COS-1, and rat primary hepatic cells (Baqué et al., 1997). Furthermore, glycogen synthase (Fernández-Novell et al., 1997) and glycogen granules (Prats et al., 2018) are associated with actin. Our results support these previous findings and suggest that the actin scaffold is the site for the initiation of glycogen synthesis. In addition, the actin scaffold was associated with endosomes, which are reportedly the location for Akt phosphorylation (Nazarewicz et al., 2011) as well as the branchpoint for Akt2-specific signaling to glycogen synthesis (Braccini et al., 2015). These findings indicate the possibility that not only actin scaffolds, but also endosomes, play roles in enhancing signals for glycogen synthesis.

Second, we also hypothesized that actin scaffold formation is involved in insulin-dependent protein synthesis, and that the process is dependent on mTOR (pSer2448) localization to actin scaffolds. Our data suggest that insulin-dependent protein synthesis may be negatively affected in CK666-treated cells through the inhibition of S6RP-mediated ribosomal biogenesis (Figure 7); although the effect of CK666 in whole protein synthesis remains to be elucidated. It is interesting that retaining mTOR to actin scaffolds would ensure not only efficient activation/phosphorylation of mTOR itself, but also effective signaling of S6RP phosphorylation. The accumulation of mTOR to actin scaffolds may enable more efficient signaling to canonical protein synthesis. Thus, aberrant translational switching centered on mTOR (pSer2448) in actin scaffolds was inferred by our PLOM-CON analysis (Figure 7F).

The actin scaffold is comprised of concave (Ω-shaped) and invaginated membrane structures, likely several fused early endosomes, surrounded by an actin mesh (Figure S1A). Liu et al. (2018) reported that activated Akt2 (pSer473) localized to phosphatidylinositol-3,4-bisphophate[PI(3,4)P2] on early endosomes, which regulates GSK3β at that location. We suppose that the concentration of activated Akt (pSer473), PI(3,4)P2, and actin-binding proteins such as SNX9 and Arp2/3 by a similar mechanism may induce the association of downstream insulin-signaling molecules, such as GSK3β (pSer9), mTOR, and mTOR (pS2344), to the actin scaffold, thereby initiating spatial-dependent regulation of glycogenesis or insulin-stimulated protein synthesis. Interestingly, PI3P on early endosomes may also be involved in this process because PI3P functions in recruiting actin-nucleating proteins, such as SNX9 (Daste et al., 2017) and WASH (Singla et al., 2019), and we also previously showed a role of PI3P on actin scaffold formation (Kano and Murata, 2019).

PLOM-CON analysis coupled with graph clustering is a powerful, new method to investigate the function of membrane-associated microdomains, including actin scaffolds as well as 2D- and 3D-condensates, for signal transduction. The power of this approach lies with its capacity in extracting a broad and varied scope of molecular information—the spatiotemporal changes in the quantity, localization, and modification states of different proteins—crucial for elucidating microdomain function.

Limitations of study

As this study was a proof of concept, a number of limitations or considerations remain to be addressed regarding PLOM-CON analysis. One such consideration is the scope of proteins that can be analyzed: The number of nodes in the network is limited on the number and characteristics of antibodies used in PLOM-CON analysis. To address this issue, we adopted two strategies. The first was to perform a pathway enrichment analysis, and second, we expanded the covariation network by combining it with a PPI database. Performing a pathway enrichment analysis in our study enabled identification of the functional pathway of “insulin resistance,” a finding which should change when actin scaffolds were disrupted. This approach, based on the transcriptome dataset, enables us to select the protein to be analyzed prior to the PLOM-CON analysis. Furthermore, the covariation network can be expanded by connecting the network with a PPI database. Using this approach, we were successful in identifying a novel protein GYG2, as a factor which acts at the actin scaffolds (Figure 5C). Thus, the implementation of both approaches enabled our PLOM-CON analysis to more comprehensively include a greater number of proteins. Another consideration is that the edge between nodes infers the correlation between the feature quantities of proteins, all of which should be experimentally evaluated. However, for several feature quantities concerning morphology (e.g., circularity and elongation), experimental evaluation may be difficult because we do not have the necessary experimental tools to strictly manipulate the morphology of intracellular structures.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Masayuki Murata (mmurata@bio.c.u-tokyo.ac.jp or murata.m.ah@m.titech.ac.jp).

Materials availability

Materials are available upon reasonable request. All unique reagents generated in this study are available from the Lead Contact with a completed Materials Transfer Agreement.

Data and code availability

The code with a toy dataset is publicly available in Mendeley Data: https://dx.doi.org/10.17632/vzz7xpymxd.1.

Experimental model and subject details

Cell culture

Rat hepatoma-derived H4IIEC3 cells were obtained from the American Type Culture Collection (ATCC), and cultured in DMEM (Nissui), supplemented with 20% horse serum (ATCC), 5% fetal bovine serum (FBS; Corning), and penicillin/streptomycin (Gibco).

Methods details

Reagents and antibodies

CK666 was obtained from Abcam (#ab141231) and 1 mM stock solution was prepared in DMSO. DAPI (1 mg/mL solution) was obtained from DOJINDO. Insulin solution (human) was obtained from Sigma–Aldrich. Rhodamine-conjugated phalloidin was purchased from Thermo Fisher Scientific. The primarily and secondary antibodies used in this study are listed in Table S1.

Transfection of LifeAct-EGFP plasmid and time-lapse imaging of transfected cells

The plasmid encoding LifeAct-EGFP was produced as described previously (Kano and Murata, 2019). H4IIEC3 cells were grown on glass-bottomed dishes (Iwaki, 4970-011) coated with collagen Type-IC. The semi-confluently cultured cells were transfected with plasmid of Lifeact-EGFP using Lipofectamine 2000(Thermo Fischer Scientific) according to the manufacturer’s instruction. After 24 h, the cells were incubated in serum-free DMEM for 16 h. After washing with PBS twice, the cells were changed to Hank’s Balanced Salt Solution, transferred on the stage of A1 confocal microscope (Nikon), and incubated for 10 minutes. The cells were stimulated by 1 μM insulin, and the time series of transfected cells’ fluorescence images were acquired at every 1 min for 60 min.

Indirect Immunofluorescence Method for PLOM-CON analysis for insulin-treated H4IIEC3 cells

H4IIEC3 cells were grown on gelatin (WAKO)-coated, 96-well, clear, flat-bottom imaging microplates (Falcon, 353219) and serum-starved for 16 h. The cells were pretreated with 50 μM CK666 or DMSO for 15 min and subsequently stimulated with 1 μM insulin in a culture medium containing 50 μM CK666 or DMSO for 1, 5, 10, 15, 20, 30, 40, 50, or 60 min at 37˚C at 5% CO2. The cells were washed twice with PBS, fixed with 3% paraformaldehyde for 20 min, permeabilized with 0.2% Triton-X 100 (Wako) in PBS for 15 min and blocked with 3% BSA (Equitech-Bio) in PBS for 30 min. Then the cells were incubated with the primary antibodies containing anti-GAPDH antibody as a cytoplasm marker for 2 h, and subsequently with the secondary antibodies containing Hoechst 33342 (DOJINDO) as a nuclear marker and rhodamine-phalloidin as an actin marker (Molecular probes).

Pathway enrichment analysis

Microarray analysis was performed on H4IIEC3 cells treated with DMSO for 75 min, 50 μM CK666 for 75 min, DMSO for 15 min followed by insulin for 60 min, and CK666 for 15 min followed by insulin for 60 min. Gene expression levels were compared between the following groups: DMSO, CK666, DMSO+insulin, and CK666+insulin. Genes whose expression levels changed >1.5-fold with a q-value of <0.1 were extracted. The q-value library of R was used to calculate the q-values. To focus on the pathways that were substantially changed by insulin stimulation in CK666-treated cells, the genes extracted from the comparison between the DMSO and CK666 treatment groups were excluded from those obtained from the comparison between the DMSO+insulin and CK666+insulin groups. We performed the over-representation test using cluster Profiler. The enriched KEGG orthology was used to extract pathways recorded in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database with a p-value <0.05 (adjusted by the BH method) and q-value < 0.1.

Automated confocal well plate imaging

The immunostained cells on 96 well imaging microplate were observed with confocal microscope (Nikon TiE inverted stand with an A1R Galvano scan head and fiber-coupled four line excitation laser LU-NV combiner equipped with 408nm-, 488nm-, 561-nm and 647nm-lines). A 40× dry objective (Plan Apo λ 40x NA0.95; Nikon), a motorized nosepiece up/down movement and a motorized stage with encode were used. 4 Channels were imaged in order of blue, yellow/orange, green and far/infrared. The gain and laser power were determined for each protein to avoid exceeding the upper fluorescence intensity limit at all time points. The hexagonal confocal pinhole was set to 1.2 airy disc size of the longest wavelength imaged. Microscope was controlled by the NIS-Elements Ver4.4 software (Nikon) to obtain Z stack images of 12 steps (1-um per step) for 9 different points per well. By obtaining 12 Z-stack images, we could cover all the cell areas in height at all time points, which ensured obtaining focused images.

With the above imaging method, DMSO image set and the CK666 image set were acquired from DMSO sample and CK666 sample subsequently with 9 points for each well, ten time points for each protein and fifty proteins in each image set.

Segmentation and 3D reconstruction

Microscopic images of the cells were analyzed using NIS-Elements Ver4.4 (Nikon) software. The nucleus was first detected from DAPI-stained images as a seed, and the cell area was determined using watershed algorithm with GAPDH-stained images. GAPDH is a housekeeping enzyme (Barber et al., 2005) and responds poorly to insulin stimulation at the protein expression level. In fact, the kinetics of GAPDH’s mean fluorescent intensity showed random patterns in all 50 time independent experiments. GAPDH’s fluorescence intensity at each time-point changed in the range of ±25% of that at 0 min, thus leading us to define the intensity’s noise range. Therefore, the feature quantities, including at least one value that changed more than or less than 25% compared with that at 0 min, were used to calculate the correlation. The cytoplasm area was then calculated by subtracting the nucleus area from the cell area. The actin domain was detected by a Detect Regional Maxima algorithm, from phalloidin-stained images, since the domain existed as bright round spots in the cell. As for Akt, Akt(pSer473), Akt(pThr308), EEA1, FoxO1(pSer256), GSK3b(pSer9), IR(pTyr1185), and WASH1, there were aggregated bright spots, and their numbers and volume were measured and included as feature quantities.

As cellular dynamics can affect the segmentation quality, the segmentation algorithm was carefully optimized and validated by examining the images for all experimental condition. Although insulin addition caused the dynamic change in actin dynamics and induced the formation of the actin domain and filopodia and lamellipodia, we were able to precisely segmentate the nucleus and cell region using the strong fluorescence signal of DAPI and GAPDH staining.

Hereinafter, the cell area for each sectioning image is referred to as a 2D cell area, and other areas (nucleus, cytoplasm, actin domain, aggregates) are referred to as 2D intracellular compartment areas. We created a 3D cell area by linking the 2D cell areas in different Z planes. More precisely, the 2D cell areas in the on-focus plane, defined as the plane where the sum fluorescence intensity of GAPDH in the cell, was maximal, and the planes above and below the on-focus plane were linked to create a 3D cell area. We confirmed that each 2D cell area in the different plane partly overlapped each other when they were projected in the same plane. To construct the 3D intracellular compartment areas, we first determined which cell area the 2D intracellular compartment area belonged to, and linked the 2D intracellular compartment areas in the same cell from the above-mentioned 3 planes to create the 3D intracellular compartment areas.

Quantification of images

In the 3D cell area and the 3D intracellular compartment area (nucleus, cytoplasm, actin domain, bright protein region), we measured the mean fluorescence intensity of each protein, the volume of each area, and the number of aggregated structures. These values, referred to as “feature quantities,” are shown in Table S2. The synthesis mean intensity is calculated by the formula, , where “” is the -th intracellular compartment in a cell “”, is the total intensity in the -th intracellular compartment, and is the area of the -th intracellular compartment. The mean intensity is calculated by the formula, , where “” is the -th intracellular compartment in a cell “”, is the mean intensity in the -th intracellular compartment, and is the number of the intracellular compartments in a cell.

Covariation network inference

Let us consider the median variables that obtained above procedures as random variables generated from the multivariate Gaussian distribution , where and represent feature and time indices, respectively. The is called the precision matrix and corresponds with the inverse covariance matrix, i.e., . If and only if an entry in , then the two variables are conditionally independent, so that there is no edge in a graph between the th and th features. Our first goal was to obtain the sparse precision matrix. The graphical lasso model can be represented as the following optimization problem:where is defined to be the element-wise L1 norm and the empirical covariance matrix is defined to be . The regularization parameter of controls the strength of the sparsity of the edges in the precision matrix. In other words, while the case that corresponds to the standard GGM, and the larger the value of the regularization parameter is, the fewer the number of edges of the network is.

We optimized the precision matrix by using the coordinate descent algorithm (Friedman et al., 2008; Ma et al., 2007; Idé, 2009; Montastier et al., 2015; Varoquaux et al., 2012). We constructed a partial correlation matrix P through a normalization to make the diagonal entries of the precision matirx with a minus, i.e., the (i,j)-th component of the partial correlation matrix was defined as

We also defined the adjancy matrix as the non-zero, off-diagonal entries in the partial correlation matrix corresponding to a graph. Note that the partial correlations can indicate the relation among all of the feature quantities in each protein. Thus, the partial correlation matrix can be represented by a block partitioned matrix per proteins and intracellular compartments. To discover the biological knowledge and to make the interpretation easier, the network was consolidated by the following procedures: 1) divide the partial correlation matrix into block matrices, 2) extract the maximal value of each block matrix, 3) reconstruct the matrix with such values, which have the dimension of the total number of the compartments of each protein (see Figure 2). In our study, we call the network depicted based on the adjacency matrix of such consolidated matrix covariation network.

In our numerical experiments, the networks were estimated by assigning several different values to the regularization parameters within 0.79 to 0.93, at a pitch of 0.02, in DMSO condition and 0.82 to 0.94, at a pitch of 0.02, in CK666 condition Then, we determined the minimum value of the under the condition that the edges between insulin receptor (IR) and FOXO1’s phosphorylated form [p-FOXO1(Ser256), p-FOXO1(Ser319), p-FOXO1(pThr24)] disappear so that the similarity of the time series of the intermediates was effectively inferred. Detailed procedures, including schematic pictures, are denoted in Figure 2.

Graph clustering

We regarded the covariation network as undirected graph G=(V,E) with n vertexes and m edges where they correspond to the compartment of each protein and the non-zero component of the matrix, respectively. We applied the Overlapping Cluster Generator (Becker et al., 2012) to the graph for decomposing a network into overlapping small networks. The OCG algorithm is the solver of an optimization problem in community detection based on an extension of Newman’s modularity function. We set initial condition as centered clique in the present work, which is known to show the high sensitivity comparing to other initial settings regardless of the network density. This is because our formulation for inferring the network can change the density of it according to the regularization parameter. For performing the method, we used R package linkcomm (Kalinka and Tomancak (2011)).

After applying the OCG, we computed the mean degree 2|E_i |/|V_i | of each cluster C_i=(V_i,E_i) and whole graph under the condition of DMSO or CK666 for each L1 regularization parameter (Figures 4A and 6A). Finally, we picked up the clusters of which the mean degrees exceed it of whole graph (Figures 4B and 6B). The subgraphs, i.e., clusters, obtained as a result of a series of the procedure are considered to be sets of functionally relevant proteins.

Validation of antibodies

Serum-starved H4IIEC3 cells were treated with or without 100 nM insulin for 15 min, and fixed, permeabilized and immunostained by antibodies against indicated proteins. Microscopic images were validated whether the appropriate staining patterns for the protein were observed. The images of the cells were analyzed using NIS-Elements Ver4.4 (Nikon) software. The antibodies with which immunofluorescent images showed sufficiently high signals (S) compared with the background (N) (i.e., where the S/N ratio was 5 or greater), were used for the analysis. We measured the mean fluorescence intensity of each protein in the areas of cell, nucleus, and cytoplasm at a single cell level, and the mean level of fluorescence intensity in the cell area and the ratio of fluorescence intensity between the nucleus and cytoplasm of the proteins were shown as a box plot.

Isolation of RNA and RT-PCR

Serum-starved H4IIEC3 cells were pretreated with DMSO or 50 μM CK666 for 15 min, and further with or without 1 μM insulin for 60 min. Total RNA was purified from the cells using an RNeasy Mini Kit (Qiagen) and reverse-transcribed using a ReverTra Ace qPCR RT Kit (Toyobo). Real-time PCR was carried out using Fast SYBR Green Master Mix (Applied Biosystems) and a StepOnePlus Real-Time PCR System (Applied Biosystems). The primer pairs for rat G6PC and PPIA are as described in Kano et al. (2017). PPIA was used as the internal standard. Three independent experiments were performed. The means and the standard errors of the relative expression level of G6PC are shown in the graph.

Western blotting

Western blotting was performed as described previously (Kano et al., 2012). Three independent experiments were performed.

Measurement of glycogen synthesis in H4IIEC3 cells by insulin stimulation

Serum-starved H4IIEC3 cells were grown on gelatin coated culture dish (ϕ140 mm), and were pretreated with DMSO or 50 μM CK666 for 15 min and then stimulated by 1 μM insulin for 120 min in culture medium containing DMSO or 50 μM CK666, respectively. Next, the cells were trypsinized, resuspended, and then sonicated, followed by measurement of glycogen synthesis with a Glycogen Assay Kit (Fluorometric) (CELL BIOLABS, MET-5023). The mean and the standard error of the amount of glycogen are shown in the graph. Three independent experiments were performed.

Expansion of the covariation network using a protein-protein interaction database

The covariation network was expanded in order to estimate drug-perturbed proteins. Graph clustering was performed from the covariation network of insulin-treated H4IIEC3 cells to obtain clusters as described above; GO enrichment analysis of the constituent proteins was performed on this cluster. Of the significantly enriched GOs, those that were descendants of the term “biological process” were extracted and to estimate the biological functions of this cluster. Proteins that are first neighbors in the protein-protein interactions that constituted the cluster and had at least one common GO term to the functions within the cluster were extracted. BioGRID was used as the database for protein-protein interactions.

Indirect immunofluorescence method for glycogen synthesis-related proteins

Serum-starved H4IIEC3 cells were stimulated with 1 μM insulin for 0 or 15 min. The cells were fixed, permeabilized, and blocked with BSA as described above. Then the cells were incubated with anti-GYG2 antibody (abcam, ab204784), anti-PPP1R3B antibody (Novus Biology, H00079660-B01P), anti-PPP1CA antibody (GeneTex, GTX105255), or anti-Myc antibody (MBL, M047-3) for 2 h, and subsequently with the secondary antibodies and rhodamine-phalloidin for 1 h. For staining of PPP1R3C and STBD1, we transfected plasmids encoding rat PPP1R3C with an N-terminal Myc tag or rat STBD1 with a C-terminal Myc tag to H4IIEC3 cells by lipofection two days before staining. Complementary DNA encoding PPP1R3C or STBD1 was amplified from cDNA of H4IIEC3 cells, using the synthetic primers 5′-GAGGAGGACCTGCTTATGAGCTGCACCAGGATGAT-3′ (forward primer for PPP1R3C) and 5′-TGTCTGGATCCCCGCTCATCGATAGGAGGCCAAGT-3′ (reverse primer for PPP1R3C), and 5′-CCCACCATGGCATCAATGGGCGCCGTCTGGTCA-3′ (forward primer for STBD1) and 5′-GATCAGCTTCTGCTCGTGAATCCCCCACCACCCAT-3′ (reverse primer for STBD1), and was subcloned into pCMV-Myc vector (Clontech).

Structured illumination super-resolution microscopy

Super-resolution imaging was performed with an N-SIM microscope system (Nikon) equipped with a SR Apo TIRF 100 × 1.49 N.A. oil immersion objective and an DU-897 X-9699 camera (Andor Technology). 3D-SIM image stacks were acquired with a Z-distance of 0.12 μm, covering the entire thickness of the actin domain (3.6 μm). Acquired images were reconstructed computationally using the reconstruction slice system from NIS-Elements software (Nikon).

Cross-Validation (CV)

Using the same cell line H4IIEC3 and 96-well plate, the cells were fixed 0, 1, 5, 10, 20, and 30 minutes after 1 μM insulin stimulation, and the cytoplasm, nucleus, and actin were stained. Acetyl-FOXO1, Akt, EEA1, Erk, FOXO1, GSK3b, NUP98, PAK, pErk, pFOXO1, pFOXO1S315, pGSK3b, PI3P, PIP2, SNX9, and WASH were stained by the indirect immunofluorescence method. For each 2D image in the acquired Z-stack image, the nucleus, cytoplasm, cell, actin domains, and proteins’ aggregates were detected using the same method as in this paper. The nuclear membrane region was calculated by further reducing the detected nuclear region and subtracting it from the original nuclear region. The 3D cell region and the 3D intracellular compartment region were constructed by linking the detected 2D cell region and 2D intracellular compartment region in the same method as described above. In the nucleus, cytoplasm, actin domain, protein aggregates, and nuclear membrane, we measured the mean each protein’s fluorescence intensity, the volume of each cellular component, and the number of aggregated structures as feature quantities.

We determined the regularization parameter ρ in the graphical lasso by applying CV to the data. The procedure of K-fold CV is as follows: First, we split the dataset into K groups. Second, we took the group as a hold-out or test data and the remaining groups as a training dataset. Third, we estimated the precision matrix by applying the graphical lasso to the training dataset with a specific regularization parameter ρ1. And then, we evaluated the likelihood with the estimated precision matrix on the test set. After retaining the score of the likelihood and discarding the precision matrix, we repeated the three steps for each unique group. Summarizing the scores, we obtained the model’s evaluation score with ρ1. In this trial, we split the dataset into five groups and performed the above procedure for a set of regularization parameters ranging from 0.185 to 0.195 in 0.01 increments and found the optimal regularization parameter 0.19 that maximized the evaluation score.

Validation of segmentation quality

To simplify the analysis, we narrowed down the condition of the experiment as follows:

proteins: IR, FoxO1_2, Akt, Akt_pSer473 and GSK3B_pSer9

treatment condition: DMSO (control) treatment

compartments: nucleus and cytosol segmentation areas

To explore the network’s sensitivity to the segmentation quality, we created two patterns of low-quality segmentations for both the nucleus and cytosol, resulting in five datasets for evaluation.

These datasets were prepared as follows: Standard data with optimal segmentation, data with nucleus segmentation 80% smaller than the standard, data with nucleus segmentation 120% bigger than the standard, data with cytosol segmentation 80% smaller than the standard, data with cytosol segmentation 120% bigger than the standard. We calculated the correlation for each dataset. Since the regularization parameter (ρ) 0.79 applied in the full dataset showed too few edges in this toy dataset, we provide data in lower ρ such as 0.05, 0.1, and 0.5 to enable discussion of the segmentation quality. The difference between the standard data and those with low-quality segmentations in correlation results was calculated as the “difference rate.” More precisely, the number of correlation pairs showing different results than the standard was counted and divided with the total number of pairs to obtain the difference rate.

Quantification and statistical analysis

Statistical analysis

Statistical analysis was performed by Tukey multiple-comparisons test. Differences were considered to be significant at ∗P< 0.05 and ∗∗P< 0.01. Data are expressed as the mean ±SEM. The number of independent experiments are shown in the legend.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Akt (clone C67E7)	Cell Signaling Technology	Cat#4692S; RRID: AB_915783
Akt(pSer473) (clone D9E)	Cell Signaling Technology	Cat#4060P; RRID:AB_2315049
Akt(pThr308) (clone C31E5E)	Cell Signaling Technology	Cat#2965S; RRID:AB_2255933
APPL1 (clone D83H4)	Cell Signaling Technology	Cat#3858; RRID: AB_2255933
Arp2 (clone FMS96)	abcam	Cat#ab49674; RRID:AB_867730
Arp3 (clone FMS338)	abcam	Cat#ab49671; RRID:AB_2257830
Bmi1 (clone EPR3745(2))	abcam	Cat#ab126783; RRID:AB_11127730
DNMT1 (clone D63A6)	Cell Signaling Technology	Cat#5032S; RRID:AB_10548197
EEA1 (clone 14)	BD	Cat#610346; RRID: AB_397830
eIF4E (clone Y448)	abcam	Cat#ab33766; RRID:AB_732125
eIF4E(pSer209) (clone EP2151Y)	abcam	Cat#ab76256; RRID:AB_1523534
Elk1 (clone E277)	abcam	Cat#ab32106; RRID:AB_732141
Erk (clone 137F5)	Cell Signaling Technology	Cat#4695S; RRID:AB_390779
Erk(pThr202_pTyr204) (clone D13.14.4E)	Cell Signaling Technology	Cat#4370A; RRID:AB_2315112
FoxO1 (clone C29H4)	Cell Signaling Technology	Cat#2880; RRID:AB_2106495
FoxO1(acAsp19) (clone D-19)	Santa Cruz	Cat#sc-49437; RRID:AB_2247212
FoxO1(pSer256)	SAB	Cat#11115; RRID:AB_895328
FoxO1(pSer319)	SAB	Cat#11136; RRID:AB_895330
FoxO1(pThr24)	Cell Signaling Technology	Cat#9464; RRID:AB_329842
GAPDH	abcam	Cat#ab83957; RRID:AB_1860478
GSK3β (clone 3D10)	abcam	Cat#ab93926; RRID:AB_10563643
GSK3β (pSer9)	abcam	Cat#ab131097; RRID:AB_11155148
GYG2	abcam	Cat#ab204784; RRID:AB_1860507
GYS2	abcam	Cat#ab83550; RRID:AB_1860507
HistoneH2A (clone D603A)	Cell Signaling Technology	Cat#12349; RRID:AB_2687875
IR (clone EPR5510(2))	abcam	Cat#ab131238; RRID:AB_11155955
IR(pTyr1185) (clone E351(2)Y)	abcam	Cat#ab62321; RRID:AB_943585
IR(pTyr999)	abcam	Cat#ab105180; RRID:AB_11150795
IRS1 (clone EP263Y)	abcam	Cat#ab40777; RRID:AB_775788
IRS1(pSer312)	abcam	Cat#ab138456
IRS1(pTyr612)	abcam	Cat#ab66153; RRID:AB_1140753
IRS2 (clone EP976Y)	abcam	Cat#ab52606; RRID:AB_880683
IRS2(pSer731)	abcam	Cat#ab3690; RRID:AB_2125790
MEK1 (clone E342)	abcam	Cat#ab32091; RRID:AB_776273
MEK1(pSer218_pSer222) (clone E237)	abcam	Cat#ab32088; RRID:AB_776275
MEK1(pSer298) (clone EPR3338)	abcam	Cat#ab96379; RRID:AB_10678572
mTOR (clone 7C10)	Cell Signaling Technology	Cat#2983; RRID:AB_2105622
mTOR(pSer2442) (clone D9C2)	Cell Signaling Technology	Cat#5536; RRID:AB_10691552
p70S6K (clone EPR13429)	abcam	Cat#ab186753
p70S6K(pSer411)	abcam	Cat#ab131459; RRID:AB_11156722
p70S6K(pSer424)	abcam	Cat#ab131436; RRID:AB_11154718
PDPK1 (clone EPR245)	abcam	Cat#ab109253; RRID:AB_10861761
PDPK1(pSer241) (clone EPR336(2))	abcam	Cat#ab109460; RRID:AB_10866450
PDPK1(pTyr9)	abcam	Cat#ab111863; RRID:AB_10867434
PI3Kp85 (clone M253)	abcam	Cat#ab86714; RRID:AB_1951326
RAC2	abcam	Cat#ab154711
RAC3 (clone EPR6679(B))	abcam	Cat#ab124943; RRID:AB_10976139
S6RP (clone 54D2)	Cell Signaling Technology	Cat#2317; RRID:AB_2238583
S6RP(pSer235_236) (clone D57.2.2E)	Cell Signaling Technology	Cat#4858; RRID:AB_916156
WASH1	abcam	Cat#ab157592
βTubullin (clone AA2)	Sigma-Aldrich	Cat#T8328; RRID:AB_1844090
Chemicals, peptides, and recombinant proteins
CK666	Abcam	Cat#ab141231
Hoechst 33342	Dojindo	Cat#H346-07951
Insulin solution	Sigma-Aldrich	Cat#I9278-5ML
Rhodamine-conjugated phalloidin	Thermo Fisher Scientific	Cat#R-415; RRID:AB_2572408
Critical commercial assays
Glycogen Assay Kit (Fluorometric)	CELL BIOLABS	Cat#MET-5023
Deposited data
Raw and analyzed data	This paper	N/A
Small toy dataset	This paper	https://doi.org/10.17632/vzz7xpymxd.1
Experimental models: Cell lines
Rat: H4-II-E-C3 cell line	ATCC	Cat#CRL-1600; RRID:CVCL_0285
Oligonucleotides
Primer: PCK1 Forward: GTGGGTCCTGGACACTGACT	Kano et al. (2017)	N/A
Primer: PCK1 Reverse: AATGCCTGACAAGACTCCA	Kano et al. (2017)	N/A
Primer: G6PC Forward: GGCAATGCTGGACCAAACAC	Kano et al. (2017)	N/A
Primer: G6PC Reverse: AAACGCTCCATGGCTTCCAC	Kano et al. (2017)	N/A
Primer: PPIA Forward: GGCAATGCTGGACCAAACAC	Kano et al. (2017)	N/A
Primer: PPIA Reverse: AAACGCTCCATGGCTTCCAC	Kano et al. (2017)	N/A
Primer: PPP1R3C Forward: GAGGAGGACCTGCTTATGAGCTGCACCAGGATGAT	This paper
Primer: PPP1R3C Reverse: TGTCTGGATCCCCGCTCATCGATAGGAGGCCAAGT	This paper
Primer: STBD1 Forward: CCCACCATGGCATCAATGGGCGCCGTCTGGTCA	This paper
Primer: STBD1 Reverse: GATCAGCTTCTGCTCGTGAATCCCCCACCACCCAT	This paper
Recombinant DNA
LifeAct-EGFP	Kano and Murata, 2019	N/A
pCMV-Myc	Clontech	PT3282-5
Software and algorithms
NIS-Elements ver 4.4	Nikon	RRID:SCR_014329
linkcomm ver 1.0-14	Kalinkaand Tomancak (2011)	N/A