GSK3β-driven phosphorylation of ABLIM1 regulates its interactions with titin cardiac muscle

The sarcomere is the main contractile element of cardiac muscle (Fig. 1). Its tens of proteins are tightly coupled to the myofilament, where they are directly responsible for generating force in response to elevated intracellular Ca²⁺. Although the structure and function of many sarcomere proteins have been elucidated, myofilament-associated proteins with intrinsic disorder (MAPIDs) (Sun and Kekenes-Huskey., 2023) are among those that remain poorly understood. Much of this challenge stems from the scant availability of biophysical data for the unstructured regions. The mechanosensitive ABLIM1 protein is an example of a MAPID serving as a key mechanosensitive regulator of muscle contractility (Anderson et al., 2021).

Actin-binding LIM protein 1, encoded by the ABLIM1 gene, is one of more than 60 proteins with LIM domains, which play crucial roles in the regulation of developmental pathways and cellular responses to mechanical stimuli (Sala and Oakes, 2023). This is done by facilitating protein–protein interactions with auxiliary proteins. Thus, the LIM domains of the ABLIM1 N terminus have the potential to mediate interactions between actin filaments and cytoplasmic targets (Fig. 2). ABLIM1 protein (88 kD) is known to traverse the actin cytoskeleton (Roof et al., 1997). Recent evidence suggests that it also establishes a connection between the Z-disk binding domains of titin, which is believed to contribute to the regulation of length-dependent activation (LDA) in cardiomyocytes (Stachowski-Doll et al., 2022).

ABLIM1 comprises four LIM domains linked to a villin headpiece by a largely intrinsically disordered region (IDR). Thus, like many other sarcomere proteins, ABLIM1 IDRs are subject to posttranslational modifications (PTMs), which suggests their importance in the rapid regulation of contraction (Sun and Kekenes-Huskey, 2023). However, their contributions to myofilament function have not been fully resolved. An indication of its role comes from a study showing that an IDR linker in a similar LIM domain–containing protein binds to plant actin to facilitate aggregation, thus influencing the remodeling of F-actin (Thomas et al., 2008). Therefore, we anticipated that IDRs in ABLIM1 may play analogous roles in regulating important actin-dependent functions in mammals.

ABLIM1 is known to be highly phosphorylated (Dong et al., 2020; Roof et al., 1997), which is confirmed by the PhosphoSitePlus database (Hornbeck et al., 2014). A study (Stachowski-Doll et al., 2022) discovered that the state of phosphorylation of ABLIM1 affects its interactions with the Z1Z2 segment of titin, which coincides with changes in myocyte contractility. Two of these sites are located between the LIM2 and LIM3 domains (see the magenta region in Fig. 3).

GSK3β is a critical kinase that phosphorylates ABLIM1, facilitating strong force generation at extended sarcomere lengths. This phenomenon is known as length-dependent activation (LDA) in healthy cardiac tissue (see Fig. 1). Its knockout in a mouse model was shown to reduce the phosphorylation of ABLIM1 at S354, S494, and S496 (see Table S1), with additional phosphorylated protein fragments remaining unchanged. The study by Stachowski-Doll et al. (2022) identified that in normal cardiac function, constitutive GSK3β activity sustains high levels of ABLIM1 phosphorylation, negligible titin/ABLIM1 interactions, and pronounced LDA. In contrast, a decrease in GSK3β activity typical of heart failure (HF) (Kirk et al., 2014; Stachowski-Doll et al., 2022) results in lower levels of phosphorylated ABLIM1, increased titin/ABLIM1 binding, and reduced LDA. Therefore, we hypothesized that the phosphorylation status of ABLIM1 impacts the ability of its LIM domains to access titin in the myofibril, where it can hinder LDA.

In this investigation, our objective was to identify how phosphorylation affects the properties of ABLIM1 in binding the Z1Z2 segment of titin. To investigate these properties, we constructed a molecular model of ABLIM1 by combining simulations of all-atom protein structures and coarse-grained (GC) representations. The all-atom simulations were performed to assess the secondary structure and key interactions that stabilize the individual LIM domains. Based on these properties, we constructed a CG model of the N-terminal domain of ABLIM1 (residues M1-S450). The CG simulations allowed us to determine the statistical properties of the ABLIM1 conformation ensemble, including the openness of the folded LIM domain regions (I97-E343) and the disordered domains (residues M1-V96 and E343-S450). Observations were interpreted with a static model of myofilament contraction. Together, these simulations provide important insights into how the structure and dynamics of ABLIM1 could influence LDA, while more broadly elucidating mechanosensitive signaling of MAPIDs as a whole.

FASTA sequences for ABLIM1 and ABLIM2 (human) were obtained from UniProt with access IDs O14639 and Q6H8Q1, respectively. The structure derived from nuclear magnetic resonance (NMR) for ABLIM2 LIM domain 2 was obtained from the Protein Databank (PDB accession no. 1V6G). The structure of ABLIM2 LIM domain 2 was derived from NMR and obtained from the Protein Databank (PDB 1V6G). A predicted structure from AlphaFold (Jumper et al., 2021) of the full-length ABLIM1 protein (human) was also downloaded from UniProt (accession ID O14639).

The ELM database was used to query the ABLIM1 sequence for potential short linear interaction motifs (SLIMs), which are sequences associated with the formation of protein–protein interactions and are often subject to PTMs. The ABLIM1 FASTA sequence (from UniProt O14639) was provided as input. The BLASTp (Altschul et al., 1990) sequence alignment algorithm was used to align all LIM domains with one another to ensure sequence identity. Each LIM domain (ABLIM2 LIM2 and ABLIM1 LIM1–4) was treated as an independent sequence to align with all other amino acid sequences of the LIM domain.

Tandem mass spec was used by Stachowski-Doll et al. (2022) to identify phosphorylated fragments from control and GSK3β KO myofilament-enriched myocardial tissue from mouse. In that study, phosphorylation was detected at ABLIM1 sites S215, S216, S354, S396, S397, S399, Y401, S407, S411, S466, S470, T473, S475, T477, S479, S494, T495, S496, S499, T586, S658, S670, S671, S752, T753, and 1858 (see Table S4). Only sites S354, S494, and S496 were found to be significantly different between samples (P ≤ 0.05), with each site having a roughly twofold reduction in the GSK3β KO tissue (Table S1).

The human sequence is 82% identical to that of the mouse, including a gap where the mouse S354 phosphorylation is found. Therefore, we modeled the equivalent positions of S450 and S452 in the human sequence (see the equivalent residues in Table S4). The protein sequences of the two species are identical in this region. Since the mouse S354 site region adjacent to LIM4 is not found in humans, we instead modeled phosphorylation of the MSSSPKET sequence found between LIM domains 2 and 3, and specifically S216, which is conserved in both species. We selected S216 because the MSSSP at site S216 (see the GSK3β label for ABLIM1 LIM2 in Fig. 3) is a common recognition motif for GSK3β and other kinases per ELM (Kumar et al., 2023). This site has been identified as a phosphorylation target in humans (PhosphoSitePlus database; Hornbeck et al., 2014; Beurel et al., 2015). In addition to phosphorylation of the region identified by Stachowski-Doll et al. (2022), phosphorylation of S215/S216 has been detected in stem cells and tumors/cancers (Klammer et al., 2012; Mertins et al., 2014), while a model of MLCK-null endothelial cells was found to reduce S216 phosphorylation (Isobe et al., 2020).

We performed atomistic molecular dynamics (MD) simulations of the isolated LIM2 domain of ABLIM1, the isolated LIM2 domain of ABLIM2, the LIM2-LIM3 domains of ABLIM1, the phosphomimetic LIM2 domain of ABLIM2 (S146D), and the phosphomimetic LIM2-LIM3 domains of ABLIM1 (S216D). The starting structures for these atomistic simulations were the resolved LIM2 domain structure of ABLIM2 (PDB 1V6G) or the AlphaFold-predicted full-length structure of ABLIM1. The missing Zn²⁺ cations in the AlphaFold-predicted structure were added by first aligning PDB 1V6G to each individual ABLIM1 LIM domain and then copying and pasting Zn²⁺ cations into the AlphaFold-predicted structure. The structure was then solvated in a TIP3P water box with a 12 Å margin and filled with 0.15 M NaCl salts. The CHARMM36m (Huang et al., 2017) force field was used to describe the system. The prepared system was first subjected to the minimization of energy for a maximum of 50,000 steps of the steepest descent algorithm (although the minimization would end early when the system reached the desired maximum force of <1.0 × 10³ kJ/mol nm²), followed by a heating stage during which the system was heated to 310 K over 2 ns in the NVT ensemble. During minimization and heating, restraints were applied to the protein backbone atoms with a force constant of 1 kJ/mol nm². To maintain coordination between Zn²⁺ and coordinating residues, harmonic restraints were introduced between each Zn²⁺ and the side chains of the coordinating residues in ABLIM1. A force constant of 1 kJ/mol nm² was applied for minimization, increased to 10 kJ/mol nm² during the NVT equilibration step, and finally reached 100 kJ/mol nm² for the final equilibration step and throughout production. A time step of 2 fs was used in most MD simulations, except for equilibration, where a 1-fs time step was necessary for simulation stability. The cutoff for nonbonded van der Waals interactions was set at 1.1 nm. Electrostatic interactions were handled using the reaction field method, with a dielectric constant of 15 within the cutoff distance and infinite beyond. The integration of Newtonian equations of motion used the velocity Verlet algorithm. The LINCS algorithm (Hess et al., 1997) was used for constraint calculations, and the velocity-rescaling thermostat and the Parrinello–Rahman barostat are used. The velocity-rescaling thermostat had a friction coefficient of 1 ps⁻¹. The Parrinello–Rahman barostat in GROMACS simulations used a compressibility constant of 3.0 × 10⁴ 1/bar. The selection of Langevin over velocity-rescaling methods was driven by two factors: (1) the preservation of the canonical ensemble and (2) ensuring that the velocity distribution maintains the correct Gaussian distribution. As for the barostat, Berendsen’s equation of motion was employed to govern the evolution of the computational box’s volume and shape. Atomistic MD simulations were performed using the GROMACS2024 package (Abraham et al., 2024). We additionally performed atomistic MD simulations on the four LIM domain construct of ABLIM1, in the wild type (WT), the phosphomimetic form (S216D), and the actual phosphorylated form (pS216), each for 500-ns production run (Sect. S1.4.2). These simulations confirmed that using phosphomimetic had a similar impact as phosphorylation on the structural properties of the four LIM domain construct, both making it more extended than the WT case. Importantly, we showed that by introducing restraints, all eight Zn²⁺ ions were well coordinated by four atoms each, reaching the theoretical coordination number of 4 for Zn²⁺ (Fig. S1), proving the effectiveness of using restraints to maintain Zn²⁺ coordination.

To enhance sampling efficiency for large ABLIM constructs, we performed CG MD with the Martini3 (Souza et al., 2021) model. Coarse graining of the atomistic structure was conducted using the martinize.py script, and the DSSP (Kabsch and Sander, 1983) program was utilized to identify the secondary structure of the atomistic model. The CG structure was solvated with a standard Martini water model using pre-equilibrated waters from the Martini website. Elastic network constraints were applied to individual LIM domains with force constants of 500 kJ/mol nm², and the lower and upper elastic bond cutoffs were set to 0.5 and 0.9 nm, respectively. For systems containing purely IDRs, no elastic network constraints were applied. Constraints were imposed between Zn²⁺ ions and the side-chain beads of the coordinating residues of the CG models following the methodology of the atomistic simulations. The system was first energy-minimized for 40,000 steps using the steepest descent algorithm and then heated to 300 K in the NVT ensemble over 1 ns. The equilibrated system was subject to a 2,000-ns production MD in the NPT ensemble at 300 K with a 20-fs time step. Temperature and pressure control were maintained using the Parrinello–Rahman V scale barostat. CG simulations were conducted utilizing the GROMACS2024 package (Abraham et al., 2024).

LDA was modeled using MATMyoSim (Campbell, 2014), a software package for simulating the mechanical properties of half-sarcomeres using cross-bridge distribution techniques, implemented in MATLAB (The MathWorks Inc., 2022). These files and parameters of the tension-pCa demo were utilized to simulate the conditions of WT and mutant myocardium in Campbell et al. (2018) and were adapted further to produce simulation data trending with the observed experimental data from Stachowski-Doll et al. (2022). The parameters modified from the default tension-pCa MATMyoSim values, but consistent across all simulation models we performed, can be found in Table S10. The parameters that were modified for each simulation model can be found in S9.

Additional supporting analyses, intermediate results, and parameter sets that support our main findings can be found in the supplemental data. Fig. S1 shows coordination of Zn²⁺ calculated from 500-ns atomistic MD of the WT four LIM domain construct. Fig. S2 shows conservation between LIM domains in ABLIM1. Fig. S3 shows predictions of disorder in ABLIM1 using the GlobPlot interface (Linding et al., 2003). Fig. S4 shows an average number of water molecules within 6.5 Å of the hinge motif (Res 212–220) in CG ABLIM1 models. Fig. S5 shows atomistic MD runs of WT and phosphorylated ABLIM1 four LIM domain constructs. Fig. S6 shows comparison of the double phosphoserine protein with the existing WT and phosphomimetic constructs. Fig. S8 shows predictions of $〈re〉$ with respect to the number of amino acids (N), based on fitting to fragment simulation data using Eq. 1. Table S1 shows phosphorylation sites from Stachowski-Doll et al. (2022) (mouse). Table S2 shows summary of parameters and distances used for the Flory polymer model. Table S3 shows table of SLIMs serving as kinase targets, as suggested in the ELM database (Kumar et al., 2024) for human ABLIM1. Table S4 shows tabulation of phosphorylated fragments identified in Stachowski-Doll et al. (2022), compared with the homologous positions in human ABLIM1. Table S5 shows parameters used in the polymer physics model. Table S6 shows data table for sites, sequences, and values of f_p and f_n. Table S7 shows bead labeling and data for phosphoserine residues in Martini3 CG simulations. Table S8 shows bonded parameters for phosphoserine residues in Martini3 CG simulations. Table S9 shows LDA parameters for MATLAB simulation of stress impacts on calcium dependent force production. Table S10 shows parameters that were changed from the default tension-pCa demo values in the model file for MATMyoSim.

ABLIM1 comprises four LIM domains (see Fig. 2) and a single villin headpiece. Each LIM domain is stabilized by two bound Zn²⁺ ions and consists of four β-strands (Fig. 2 C). The LIM domains exhibit 20–30% sequence identity (see Figs. 3 and S2), with highly conserved cysteine and histidine residues that coordinate Zn²⁺. Conserved glycines may confer more conformational variability to the LIM domains, given their less restrictive ϕ and ψ angles.

Several tools predict that many regions of the ABLIM1 sequence are intrinsically disordered. To illustrate this disorder, we show predictions from IUPred2A (Mészáros et al., 2018) in Fig. 2 A and similar results from GlobPlot (Linding et al., 2003) in Fig. S3. This analysis reveals that almost 30% of the sequence is disordered, primarily within its N terminus (M1-V96), as well as the segment spanning LIM4 and the villin headpiece (Q338-M696). Exceptions include folded regions within I97-E343, where its four LIM domains are located, and the C terminus beginning with position 710, where its villin headpiece is found.

To gain insight into the topology of IDRs, we used the qualitative phase diagram from the Pappu laboratory (Holehouse et al., 2017), which relates the positive and negative fractional charge (f_p/f_n) to a polymer’s topology, like swollen coils. We first used PONDR-VLXT (Li et al., 1999; Romero et al., 2001) to predict IDR fragments in the ABLIM1 protein. Our analysis revealed several predicted IDRs distributed in the ABLIM1 sequence, with most concentrated in the C-terminal half of the protein. This IDR distribution is consistent with the disordered scores predicted by IUPred2A as shown in Fig. 2 A.

To understand how the conformational states of these IDRs would change after phosphorylation, we consulted the PhosphoSitePlus database (Hornbeck et al., 2014) to collect all reported phosphorylation sites that fall within the predicted IDRs. As shown in Fig. 4 B and Fig. 5, the IDRs of ABLIM1 are enriched with phosphorylation sites. Upon plotting these IDRs into the Pappu diagram using their fraction of positive/negative charge fraction as coordinates, the IDRs’ projections are largely concentrated around the cusp between conformational phases, namely, around the molten globules and chimeras. We then updated the f_n values based on the number of phosphorylation sites we identified from PhosphoSitePlus, which shifted nearly half of the IDRs upward along the y axis. Hence, based on the change in fractional charge, phosphorylation may shift several of the IDRs toward more extended configurations and potentially even swollen coil conformations (Fig. 4 C).

Phosphorylation and other PTMs that alter sequence charges can shift conformation ensembles between different topologies (Hornbeck et al., 2014). For ABLIM, this potentially modifies its structure or interaction partners (Sun and Kekenes-Huskey, 2023). We recomputed the fractional charge density using equation f = N_p,n/N_total, where N_p,n is the number of R/K or D/E residues, to account for phosphomimetics at three sites: S216, a putative site based on high-throughput hits (Stachowski-Doll et al., 2022), and two sites S450 and S452 in the disordered C-terminal region (see Fig. 4 B).

These shifted sequences along the negatively charged axis (see red dots in Fig. 4 C). For the three unphosphorylated sites, the IDRs containing them are located near the transition zone between molten globules and chimeras. The phosphorylation of all three sites is expected to result in the IDRs adopting more molten conformations (Fig. 4 C). We anticipate that these changes will modify the orientations of the LIM domains relative to each other and will influence the effective length of the linker between the LIM domains and the villin headpiece.

The sequence analyses in the previous section suggest that ABLIM1 exhibits significant intrinsic disorder, prompting us to perform MD simulations to characterize the structure and dynamics. To begin, we utilized the atomic resolution structure available for the LIM2 domain of ABLIM2 (PDB 1V6G) and the AlphaFold-predicted structure for ABLIM1, as ABLIM1 does not yet have a crystallographically determined structure available. The LIM2 domains of ABLIM1 and ABLIM2 share 67% sequence identity (85% sequence positives) according to a BLAST protein multiple sequence alignment (Altschul et al., 1990) of the LIM domains (see Fig. 3). Therefore, we anticipated that the ABLIM2 structure should be a reliable surrogate for ABLIM1, given their high sequence similarity. A comparison of the structures is shown in Fig. 2.

We then performed all-atom MD simulations using ABLIM2 LIM domain 2 (Fig. 3). The binding of Zn²⁺ ions, which are important for stabilizing the LIM domains (Michelsen et al., 1993), was examined. For ABLIM2, Zn²⁺ ions are stabilized through cysteine, histidine, and aspartic acid residues positioned via β sheet folds (zinc finger 1: [C18, C21, H38, C41], zinc finger 2: [C44, C47, C66, C69]). Consequently, we simulated the ABLIM1 fragment comprising LIM domains 2 and 3. The resulting data (Fig. 5) suggest that these interactions are maintained throughout the simulation by visual inspection for ABLIM1 as well, confirming the reliability of the atomistic model; namely, the distances between the bonds to zinc are ∼2.5 Å and are moderately higher than the value found for the crystal structure of LIM2 (2.3 Å).

We report the corresponding root mean square fluctuations (RMSF) in Fig. 6 A for the ABLIM1 and ABLIM2 structures to assess the flexibility of the ABLIM structures. The smallest RMSF values are around the Zn²⁺-binding sites, consistent with the protein’s contacts with zinc ions and β sheets remaining stable throughout the simulation. This suggests that the folded structure of the protein is relatively immobile in this region. The R_g (Fig. 6 B) is centered about 17 Å, demonstrating a compact and well-folded structure. Calculations of R_g for ABLIM1 include LIM domains 2 and 3; thus, the sizes are larger compared with the isolated LIM2 of ABLIM2. We note that the distances between the two Zn²⁺ ions assume values close to 22 Å, which aligns two of the β sheets within the core of the protein. This configuration may be achieved through a weak hydrogen bond (∼4.5 Å between residues D92 and K132 seen in ABLIM2 LIM2; see Fig. 6, D and E). The relative openness of this region between the two half LIM domains may indicate that the two sections could move independently of each other to bind targets. This is supported by our observation that the distribution of Zn²⁺ distance values spans 17–19 Å, suggesting that the protein is relatively mobile on the nanosecond timescales we simulated.

We next performed simulations of the S216D phosphomimetic to approximate the structure and dynamics of ABLIM1 in its physiological, phosphorylated states relative to the (unphosphorylated) WT version. We found that the phosphomimetic at S216D did not significantly impact the RMSF trends relative to the WT (Fig. 6 A), whereas the R_g decreased somewhat for ABLIM1 (Fig. 6 B). With respect to the distances between the Zn²⁺ atoms, phosphomimetics tended to shift them toward larger values, with smaller changes observed for the ABLIM1 construct relative to ABLIM2 (Fig. 6 C). This increase in distance suggests an opening of the half LIM domains in relation to the WT.

In several reports (Deane et al., 2004; Sun et al., 2020; Winkelman et al., 2020), it has been observed that more than two LIM domains are generally required for LIM domain–containing proteins to bind strained actin. As an exception, testin requires a single LIM domain to bind strained actin in cells (Sala and Oakes, 2021). Other protein–protein interactions involving LIM domain proteins are typically mediated by their LIM domains or by additional functional domains, such as PDZ domains (Kadrmas and Beckerle, 2004; Sala and Oakes, 2023). We therefore expanded our ABLIM1 model to include all four LIM domains connected by unstructured regions, in the event that multiple domains could be used for binding to the Z1Z2 segment. Since these were unstructured, we anticipated that the LIM domains could have considerable rotational flexibility that might be insufficiently sampled by conventional atomistic simulations. Therefore, we resorted to Martini CG models to enhance conformational sampling by using larger integration time steps. Thus, we created Martini CG models of the ABLIM1 LIM domains from its atomistic equivalent (see Fig. 7 A). Constraints were placed between the Zn²⁺ ions and the residues with which they coordinate; secondary structure assignments were made on the basis of results from the atomistic simulations. After equilibration and 2,500 ns of simulation, the CG representations of the LIM domains remained similar to those of the atomistic data (see Fig. 7).

We report in Fig. 7 B the structure of the CG representation’s cation binding domain, which shows coordination similar to that in the atomistic equivalents. We report these distances in Fig. 7 C, but this time using the Cα atom of each side chain of the atomistic reference for comparison. The CG model of ABLIM maintains coordination distances similar to those in the atomistic equivalents. The beads of the CG models encompass multiple atoms, so their centers do not exactly match the side-chain distances reported in the atomistic models, which explains the greater distances in CG structures. Additionally, the CG models show weakening of the secondary structure alignment, which slightly increases the distances between the coordinating residue side-chain beads and the zinc ions they coordinated. In addition, the RMSF shown in Fig. 7 D shows lower values near the folded LIM domains and especially near the location of the zinc finger residues, similar to the atomistic data. This suggests that our model is a faithful representation of the atomistic reference with regard to cation coordination.

Next, we examined the R_g of the WT and phosphomimetic Martini models (see Fig. 7 E). R_g for CG is considerably smaller than the atomistic simulations, which we attribute to the different masses and positions of atomistic Cα versus backbone beads, since the structures appear similar by visual inspection (see also panel B structures). We next describe the Zn²⁺-Zn²⁺ distances (Fig. 7 F), for which ABLIM2 LIM2 exhibited moderately smaller Zn²⁺-Zn²⁺ distances (∼15 Å) than WT (∼17 Å). In comparison, their distances in LIM2 from ABLIM1 are a few Å larger than those of ABLIM2. The ABLIM1 LIM2 region (WT_ABLIM1_156-218, shown in Fig. 7 F in pale blue) showed the most similar LIM2 Zn²⁺-Zn²⁺ distances to ABLIM2 LIM2, indicating that the presence of other LIM domains may impact the openness of a LIM domain. Furthermore, the phosphomimetics appear to increase this distance in all models compared with their respective WT counterparts, with the modified extended model that includes the IDRs on either end of the LIM domains (M1-V96 and E343-S450) showing the greatest increase in LIM2 Zn²⁺-Zn²⁺ distance. Despite the differences between the properties of the atomistic and CG structures, an overlay of the domains shows considerable similarity (Fig. 7).

We extended the four LIM domain construct to include flanking regions of roughly 100 amino acids on either side. This was done to determine whether the presence of disordered flanking regions would impact the structure of ABLIM1. Consistent with the atomistic data (Fig. 6), phosphorylation did not significantly impact the values of the RMSF, except for a moderate reduction in the lability of the disordered N terminus of ABLIM1 (see Fig. 8 A). The R_g of the WT construct was ∼25 Å, which shifted to around 28 Å upon introduction of the S216D phosphomimetic. For comparison, the WT disordered N- and C-terminal regions (M1-V96 and E343-S450) assumed R_g of ∼24 Å and 16 Å, respectively, versus the region containing LIM domains that spanned I97-E343 that was centered around 23 Å. To demonstrate the differences in R_g, we provide snapshots of the WT and S216D structures in Fig. 8 C, which show that the WT forms a more compact configuration. Hence, the phosphomimetic appeared to “straighten” the alignment of the LIM3 and LIM4 domains, leading to distances of roughly 50 Å from S216 to site E343 for the S216D construct, versus 40 Å for the WT (Fig. 8 D). We anticipated that this straightening may play a role in the binding of ABLIM1 to other target proteins. This straightening is paralleled by an increase in hydration (Fig. S4).

We next constructed additional models with two phosphomimetics (S450D and S452D) in the disordered C-terminal region to assess whether phosphorylation consistently impacts R_g for longer ABLIM1 sequences compared with the LIM1-4 construct. Simulations were then performed with the double S/T-to-D phosphomimetic variants. The results revealed that the double phosphomimetic impacted the R_g of the isolated C-terminal disordered simulations, apparently narrowing the range of conformations assumed by ABLIM1, as shown in Fig. 9 B. In addition, local differences were observed in the system, indicating that the double phosphomimetic exhibited less rigidity in the structure and allowed collapse of the structure (Fig. 9 A).

Together, these data demonstrate that the WT remains extended, whereas the N-terminal region of the IDR, housing the phosphomimetics, starts to contract. Ultimately, the phosphomimetics decrease the end-to-end distance of the IDR, as shown in Fig. 9 C. Evidently, the collapse of the chain near the double phosphorylation sites restricts the range of conformations that the IDR can adopt, affecting its end-to-end distance. Hence, these data suggest opposing effects of phosphorylation on different chain sizes, with increases generally observed for the LIM domain–containing region with S216D, versus a decrease observed for phosphorylation of sites S450 and S452 within the C-terminal IDR.

To confirm that our observed trends with phosphomimetics were consistent with a phosphorylated configuration, we performed simulations using explicitly phosphorylated serines. Atomistic and CG MD simulations, employing pSerine parameters adapted from Amber ff19SB (Raguette et al., 2024) and from Pluhackova et al.’s work (Pluhackova et al., 2021), respectively, both with a charge of −2e, showed that pS216 produces an extended conformation in 500-ns atomistic MD simulations on the four LIM domain construct (Fig. S5), consistent with our phosphomimetic findings. Furthermore, our CG simulations that evaluate the E343-S473 fragment indicate a similar reduction in end-to-end distances of stable conformations for the protein phosphorylated at pS450/pS452 compared with the WT, mirroring the phosphomimetic results (Fig. S6). These results support the use of phosphomimetics as surrogates for constitutively phosphorylated proteins (Oberoi et al., 2016).

To understand the physical basis of opposing trends, we evaluated the fractional charge density of amino acid sequences near phosphorylation sites, as previously done in Fig. 4 C. For site S216, the phosphorylatable residue is flanked by the sequence motif AQPMSSSPKETTF. We selected this number of residues given that the Debye length, within which electrostatic screening effects are minimal, is ∼1 nm at typical cellular ionic strengths (Tadmor et al., 2002). The corresponding fractions of positively (f_p) and negatively charged (f_n) residues in the WT are 0.08 and 0.08, respectively (Table 1), which counterbalance (e.g., net charge density |f_p−f_n| = 0). With S216D, f_n increases to 0.15, reflecting an overall increase in the net charge density. When there is a net charge due to unequal f_p and f_n, disordered peptides tend to have extended conformations induced by self-repulsion (Mao et al., 2010), hence explaining the straightening of the LIM domains we observed earlier (see also Fig. 4). This phenomenon is also observed for charged brush polymers that tend to swell at low ionic strength versus collapse at high salt when electrostatic screening is greatest (Shefer et al., 1993).

We applied a similar approach to examine the set of amino acids at sites S450 and S452, selecting a sequence that incorporates both residues, RMIHRSTSQGSIN. The WT form in this case had a net charge density, f_p − f_n, of 0.15, which was expected to exhibit self-repulsion. By mutating the two serines into aspartic acids, the charge density was neutralized. Hence, we expected that the neutralized sequence would contract somewhat relative to the WT, which was consistent with our simulation data for the C-terminal fragment (see Fig. 9 C), for which the diminished electrostatic repulsion compacted the IDR. Jin and Gräter (2021) also reported a similar trend for Ash1, which is a positively charged 81 residue IDR with f_p − f_n = 0.185; namely, they reported a more compact ensemble after multisite phosphorylation. They additionally compared the effects of the IDR with a phosphomimetic (−1e) versus a protonated, phosphorylated (−1e) state and a deprotonated, phosphorylated state (−2e); all cases rendered Ash1 more compact. Thus, the introduction of negative charges, be it from phosphorylation or phosphomimetics, to a positively charged IDR tends to compact its conformational ensemble.

Several lines of evidence suggest that ABLIM1 binds to titin (Stachowski-Doll et al., 2022). Coimmunoprecipitation data from Stachowski-Doll et al. (2022) support binding in vivo, while in vitro experiments further confirmed that binding occurs between ABLIM1 and titin’s N-terminal Z1Z2 segment. Moreover, they also showed that the binding is phosphorylation-dependent, given that GSK3β activity attenuates ABLIM1-Z1Z2 segment binding (Stachowski-Doll et al., 2022). Based on these findings, our simulations have so far predicted that phosphorylation (S216D) opens the four LIM domain constructs relative to the WT case using CG simulations. We confirmed this observation using atomistic simulations of an ABLIM1 construct with four LIM domains (see Supplemental material).

We therefore examined the impact of phosphorylation on ABLIM1/Z1Z2 segment binding, via rigid-body docking. To this end, we selected representative conformations of ABLIM1 for both WT and S216D cases and performed ClusPro docking of these conformations to the Z1Z2 segment. This strategy was partially justified by the fact that the Z1Z2 segment is known to exist in an extended conformation in solution (Marino et al., 2006). For each ABLIM1 conformation, we considered all poses generated by docking, of which in Fig. 10 are binding poses for the WT S1 state (the compact ABLIM1) and S216D S3 state (the extended ABLIM1 conformations). Interestingly, ABLIM1 in both cases tended to bind to the N terminus of the Z1Z2 segment, while S216 was located far from the binding interface in all poses (Fig. S7), suggesting that phosphorylation does not likely affect binding via directly modulating the interface. We speculate that such binding poses are reasonable since the N terminus of intact titin is a target for protein–protein interactions (Zou et al., 2006). Among the docking poses identified, we identified several for which a β-strand from ABLIM1 interacts with another from Z1. The β-strand stacking is similar to the β-sheet formed at the interface in the Z1Z2 segment with telethonin (TCAP) when bound as a complex (PDB accession no. 1YA5 [Zou et al., 2006]; see Fig. 11 B). We further performed MMGBSA calculations on these ABLIM1-Z1Z2 segment binding poses to evaluate their stability, which indicated that the poses exhibited favorable binding free energies (∆G ≈ −40 kcal/mol; see Fig. 10), for all poses and in a phosphorylation-independent manner.

Our simulations and analyses of WT and phosphomimetic ABLIM1 structures suggest that the charge introduced by the phosphate groups affects the compactness of IDR. Therefore, we sought to investigate how the complete ABLIM1 could interact with titin within an intact myofibril. To this end, we employed a Flory scaling model for polymers (Eq. 1) to predict the geometric properties of the IDR bridging LIM4 and the villin headpiece, namely, the R_g and the end-to-end distance. Since simulating the complete structure is not feasible, we resorted to simulations of smaller IDR ABLIM1 constructs with and without phosphomimetics to predict their R_g as a function of length. This approach allowed us to estimate the scale factor ν, with which we can scale the R_g to the full length of the ABLIM1 C-terminal IDR.

We assumed that villin binds to a region on actin directly adjacent to the terminus of an actin from a neighboring sarcomere (see Fig. 1). This binding location is consistent with the idea that villin binds to actin (George et al., 2020). The tethering of the C terminus to actin would leave the N-terminal portion of ABLIM1 free to rotate and bind to the Z1Z2 segment dimer formed by paired titins. We next compared our sequence properties with those in Hofmann et al. (2012), where they reported data on the relationship between the R_g and the number of segments (amino acids) for folded and unfolded proteins. In Fig. 12 A, we compiled R_g data from Hofmann et al. (2012) and our own simulations. Their study fit a form of Eq. 1 to describe this relationship for folded (orange) and unfolded proteins (blue), where the latter exhibited a larger R_g. The fits for the folded and unfolded proteins were ν = 0.34 and ν = 0.58, respectively (Hofmann et al., 2012).

To compare our results with the disordered region of ABLIM1, we extracted R_g estimates from our MD simulation data. For the WT fragment spanning residues 343–473, the average radius of gyration (⁠$⟨Rg⟩$⁠) of 53 Å for the 131 amino acid sequence (see Fig. 9) corresponds to ν = 0.68 in Eq. 1. Although this value deviates from the 36 Å estimate based on the Hofmann et al. (2012) fit for unfolded proteins, they also reported a larger-than-expected R_g = 30 Å for a small intrinsically disordered protein ProTα. In that case, ν = 0.65 was obtained, which more closely matches our result. For our molecular simulations of the phosphomimetic fragment that yielded a $⟨Rg⟩$ of 48 Å, we obtained a moderately smaller scale factor of ν = 0.65. This decrease is attributed to the two S-to-D mutations neutralizing the positive charge density of the WT, mitigating swelling due to self-repulsion. Using our estimates for ν based on the simulated fragments, we predicted the corresponding $⟨Rg⟩$ values for the full (366 amino acid) C-terminal sequence as 93 and 82 Å for WT and phosphomimetic forms, respectively.

To relate the geometries of the ABLIM1 fragments to their ability to interact with the Z1Z2 segment on neighboring actin filaments, we use the probability distribution P(r) from Eq. 2. This relationship determines the likelihood that a fragment assumes a distance comparable to the spacing of the actin filaments. We first compute the average end-to-end distance $〈re〉$ for the intact C terminus of ABLIM1 using predictions from our molecular simulations of shorter fragments. For this, we report in Fig. S8 the $〈re〉$ values for the short WT and phosphorylated ABLIM1 C-terminal fragments, which were found to be 170 and 132 Å, respectively, from Fig. 9 C. Using Eq. 1 again, but this time fixing ν to the values determined from the R_g fits, we solved for R₀ using the $〈re〉$ values from our simulation data. Our fit with the fixed ν values yielded R₀ = 6.0 Å for linear scaling of the end-to-end distance versus the number of residues in Fig. S8. The fitted parameters provided reasonable estimates of the $〈re〉$ for our simulated fragments, though they under- and overestimated the values to an acceptable degree. Extrapolated fits to the Flory equation Eq. 1 for the 366 amino acid C terminus resulted in $〈re〉$ values of 311 and 276 Å, respectively.

We then evaluated how the average $〈re〉$ influenced the chain distributions of the WT fragment and its phosphomimetic counterparts. These analyses yielded P(r) for each end-to-end distance, as shown in Fig. 12. To obtain these results, we substituted the estimated average end-to-end distance (⁠$〈re〉$ = 170 Å) for the WT fragment into Eq. 2, yielding a maximum at approximately r_e = 120 Å. Based on $〈re〉$ values of 340 and 230 Å for the full-length C-terminal cases, we observed that P(r) maxima were right-shifted toward larger lengths of roughly R_e = 190 Å for the phosphomimetic case versus about 270 Å for the WT.

Lastly, we compared the $⟨re⟩$ of the ABLIM1 protein when bound by its villin domain to an actin filament, relative to the spacing between the actin filaments (see the schematic in Fig. 1); namely, we estimated the end-to-end distance from ABLIM’s bound villin domain to its LIM domains to assess the LIM domains’ likelihood of interacting with an adjacent actin filament capped with titin Z1Z2 segments. The actin lattice spacing in the Z band was estimated to be ∼30.4 nm (Irving et al. 2011), which is marked in Fig. 12 B by a blue circle. The P(r) data in Fig. 12 B suggest a marginally higher probability for the WT C-terminal region versus the phosphorylated condition will assume a value R_e matches the filament spacing. This is indicated by the shaded region in Fig. 12 B. At this distance, the WT ABLIM1 is more likely to interact with a Z1Z2 segment dimer on an adjacent actin relative to its phosphorylated counterpart; by extension, the WT (unphosphorylated) case may have a greater potential to interfere with the Z1Z2 segments that anchor titin to actin. Next, we assumed that ABLIM1 is no longer perpendicular to actin but rather anchored at an angle similar to α-actinin (see Fig. 1); in this scenario, larger R_e values would be needed to bridge adjacent filaments, which results in the WT being increasingly favored over the phosphorylated ABLIM, as illustrated by the shaded region.

To reflect disrupted passive stiffness attributed to titin similar to the reduced activity of GSK3β, we reduced k_p and L₀ for our GSK3β knockout model configuration to match the experimental trends reported by Stachowski-Doll et al. (2022) (Table S9). For this configuration, we observed a comparable maximum stress for the resting sarcomere lengths, consistent with data from Stachowski-Doll et al. (2022), when coupled with a change in the k₁ rate constant from 1 to 2, although our results show a loss of cooperativity evidenced by the smaller slope of the stress/pCa curve. The similar level of stress despite a weaker passive force term is due to our assumption that the resting sarcomere length was shortened. We based our assumption on a similar treatment for titin mutants in Campbell et al. (2018). We found that simulations of the GSK3β knockout predicted a smaller increase in maximal force at the extended sarcomere length compared with WT. Consequently, the values of ∆ pCa50 were reduced by two thirds to 0.05. In other words, reducing the passive tension in our model was sufficient to reflect the loss of LDA reported by Stachowski-Doll et al. (2022), despite some quantitative differences in cooperativity, which could be attributed to the simplicity of the kinetic scheme used.

We note that the Campbell model (Campbell et al., 2018) exhibits elevated force at low pCa (<6.2) for stretched tissue, which is consistent with force/pCa data collected by Patel et al. (2012) using rat skinned trabeculae. However, the data collected by Stachowski-Doll et al. (2022) used isolated skinned mouse cardiac myocytes, for which the forces at low pCa converged to nearly zero at low pCa. Because we limited our changes to the Campbell model to the parameters governing the myosin off-to-on transition and slack length, our predicted tension for stretched sarcomeres in trabeculae at low Ca²⁺ was elevated relative to the slack condition, in contrast to the data from Stachowski-Doll et al. (2022) for isolated myocytes. We did not further refine the model given the many differences in passive tension and force development in intact tissue compared with isolated cells that would be difficult to reformulate. Nonetheless, both preparations exhibit loss of LDA when the ability to generate passive force is reduced relative to WT, thus validating our hypothesis.

In this computational study, we examined the myofilament protein ABLIM1 and its interactions with titin at the Z1Z2 segment. Our findings suggest that ABLIM1 LIM domains are compatible with binding titin and that phosphorylation can influence the conformation distribution of ABLIM1 relative to the WT (unmodified) protein. Furthermore, intact ABLIM1 is capable of spanning typical actin–actin spacings at the z-disk, and its span is sensitive to the phosphorylation state. These molecular findings provide a mechanism linking phosphorylation-dependent ABLIM1 conformation and titin interaction to the modulation of passive tension and, consequently, LDA. Overall, these data have significant implications for mechanosignal transduction in cardiac myocytes, particularly in the context of failing hearts where reduced ABLIM1 phosphorylation has been observed (Stachowski-Doll et al., 2022). These results also strengthen the support for LIM domains in mediating diverse mechanosensitive functions (Anderson et al., 2021).

We elucidated the binding partners of ABLIM1 by analyzing the STRING database. Our search of the database revealed that ABLIM1 binds to other sarcomere proteins, including titin, β-actin, and LIM-nebulette, with varying degrees of confidence (Szklarczyk et al., 2023). This finding is consistent with previous research that indicates that ABLIM proteins are generally associated with binding actin, often through their villin-type “headpiece” domain common to tens of cytoskeletal proteins (Vardar et al., 2002). Additionally, LIM domains, such as those present in ABLIM1, have been suggested to act as protein-binding modules (Dawid et al., 1998) and may recognize different actin conformations, particularly when actin is stressed. Recent MD simulations have provided a mechanism by which LIM domains can bind strained actin conformations, including the recognition of two binding sites exposed at the cracked interface (Zsolnay et al., 2024).

Intriguingly, we identified potential associations of ABLIM1 with KCNJ12 and CLCN1 by STRING analysis. KCNJ12 codes for an ATP-dependent inward rectifying potassium channel (Kir2.2), and CLCN1 is a gene for a chloride channel. ABLIM1 has previously been identified as an actin-binding LIM protein associated with Kir2.2 in brain proteins (Leonoudakis et al., 2004). Therefore, the expression of Kir2.2 in the heart (Reilly and Eckhardt, 2021) could exhibit similar interactions. A previous study indicated that CLCN1 and ABLIM1 are spliced by the same protein MBNL1 (Xing et al., 2021), but evidence for coexpression and association is lacking. The potential association between ABLIM1 and these ion channels may imply a role of ABLIM1 in modulating cell excitability. This could provide a link between mechanics and excitation at the z-disk, although more research is needed to confirm this hypothesis.

ABLIM1’s binding of titin through its N-terminal Z1Z2 segments has been observed (Stachowski-Doll et al., 2022), although the specific region responsible for this interaction remains experimentally undetermined. Our data suggest that LIM domains have the capacity to bind the Z1Z2 segment directly, as individual LIM domains can bind to the Z1Z2 segment of titin, as shown via our docking simulations and MMGBSA calculations. In fact, some of the predicted binding poses align a LIM domain β-strand with those of the Z1Z2 segment, similar to the resolved Z1Z2–TCAP complex (PDB accession no. 1YA5) where intermolecular β-strands interact tightly as antiparallel sheets (Zou et al., 2006) (Fig. 11). This binding mode is also reported for LMO4 and LDB1 (1RUT) (Deane et al., 2004), where β strands of the LIM domains align antiparallel with those of the target protein. The TCAP–Z1Z2 segment complex suggests that TCAP dimerizes Z1Z2 segments from two titin strands, bridging the β-sheets between a Z1Z2 pair (Fig. 11). In this capacity, ABLIM1 could compete with or compromise the binding of TCAP to the Z1Z2 segment. Consequently, TCAP mutations have already been established in hypertrophic and dilated cardiomyopathies (Hayashi et al., 2004). This raises the possibility that the disruption of the binding of TCAP to titin by ABLIM1 could contribute to cardiac dysfunction.

We show that phosphorylation of ABLIM1 may play a crucial role in determining its conformational ensemble. Specifically, we speculated that phosphorylation primarily affects LIM availability to the Z1Z2 segment by controlling the extension of the linker that connects the LIM domain and the villin headpiece. We propose that this sets the local concentration of the LIM domains near the titin N terminus (Fig. 11 A). Hence, this proposed model suggests that the association of the villin headpiece with actin could allow the LIM domains to bind to the Z1Z2 segment, with the potential to compete with or otherwise disrupt the interactions between TCAP and the Z1Z2 segment. Our docking results surprisingly indicated that phosphorylation between LIM domains does not have a significant impact on the interactions of the Z1Z2 segment, apparently in contrast to previous findings from an in vitro binding assay involving the full-length ABLIM protein (Stachowski-Doll et al., 2022). However, computational limitations precluded us from docking larger constructs, for which multisite phosphorylation could, in principle, disrupt its interaction with the small titin segment. It is also possible that the association of the Z1Z2 segment with ABLIM1 may also be controlled by the specific phosphorylation pattern within the ABLIM1 disordered C terminus. The ELM database suggests that ABLIM1 has numerous potential GSK3β phosphorylation sites, of which several were found to be phosphorylated by mass spectrometry (Stachowski-Doll et al., 2022). In addition, a study on the recruitment of ABLIM1 to cortical actin found that phosphomimetics at S214 and S431 were sufficient for ABLIM1 to colocalize with the cell membrane (Dong et al., 2020), evidencing that sites both within and outside the region containing the LIM domains have functional phenotypes. Changes in ABLIM1 phosphorylation could also impact LIM/actin interactions among other targets of the LIM domain and thus influence localization in vivo. Therefore, it is important to quantify the full range of potential phosphorylation sites in human ABLIM1 to appreciate the phosphorylation and GSK3β dependence of titin association.

Based on our ABLIM1 simulation data, we aimed to relate the interfilament spacing in the sarcomere to the chains’ end-to-end distance. Our calculations revealed that the LIM domains of the full-length protein would likely be found ∼82 Å from the C-terminal villin headpiece when phosphorylated versus 93 Å when not. Hence, the ABLIM1 LIM–Z1Z2 segment interaction would be most likely when the disordered domain spanning the N-terminal LIMs and C-terminal villin is not phosphorylated. In principle, if ABLIM1 were tethered to actin by the villin headpiece, the LIM domains could interact with the titin Z1Z2 segments on an adjacent thin filament, given interfilament distances of ∼30 nm reported in Irving et al. (2011). Moreover, as the sarcomere lattice contracts under isovolumetric stretching (Minozzo et al., 2013), the likelihood of interaction would further increase. This mechanism is analogous to the role of phosphorylation in toggling interactions between the N-terminal extension of MyBPC and actin, where the protein is tethered to the thick filament by its C-terminal domains (Rahmanseresht et al., 2021).

Phosphorylation’s ability to modulate interactions with the Z1Z2 segment could potentially influence myofilament mechanics. We believe that the most immediate impact would be reflected in the LDA of the sarcomere, which describes how more force is generated as the sarcomere is stretched. This phenomenon is a fundamental part of the Frank–Starling effect, where fewer equivalents of calcium are needed to achieve equivalent numbers of cross-bridges relative to when the sarcomere is at its resting length. A recognized contribution to LDA is the passive tension arising from titin, which helps restore the filaments to their equilibrium length when stretched (Campbell et al., 2018). It has been suggested that this passive tension enhances the recruitment of myosin heads in the “off state” to the “on state”; the myosin on state facilitates the formation of force-generating cross-bridges with the thin filament (Campbell et al., 2018). We postulate that some of this passive tension is due to titin binding to actin termini through dimerization of their Z1Z2 segments through TCAP, which is supported by Pinotsis et al. (2006). In turn, factors that attenuate this anchoring are expected to decrease LDA, as has already been argued for titin treatments that disrupt passive tension (Tombe et al., 2010). Mutations in the Z1Z2 segment have also been shown to reduce its affinity for TCAP by up to 40%, which would likely impair passive tension (Itoh-Satoh et al., 2002). Hence, we speculate that interactions (such as the association of the LIM domain) that interfere with the binding of the Z1Z2 segment to TCAP could also have an impact on LDA. Along these lines, TCAP deficiency leads to animals being less adaptive to stress due to reduced fractional shortening (Knöll et al., 2011).

Using the Campbell myofilament contraction model, we fit its parameters to resemble the pCa50 data from Stachowski-Doll et al. (2022). This was done in part by reducing the myosin off-to-on transition rate owing to reduced passive tension for the GSK3β KO case. Our simulation results showed a markedly reduced ∆pCa50 in the GSK3β knockout condition (0.05 compared with 0.26 in WT). The reduction in ∆pCa50 supports the hypothesis that decreased passive tension underlies the observed loss of LDA. Although the study by Stachowski-Doll et al. (2022) did not explicitly examine the effects of ABLIM1 phosphomimetics on LDA, evidence from MYBPC phosphomimetics and troponin I suggests that aspartic acid substitutions could reproduce the effects of phosphorylation on LDA (Kumar et al., 2015).

A direct connection between ABLIM1 regulation and cardiac dysfunction remains to be established. However, in a canine model of HF, decreased GSK3β activity was reported, which was associated with a reduction in ABLIM1 phosphorylation (Kirk et al., 2014). Similarly, ABLIM1 knockdown in C2C12 myoblasts has been shown to decrease serum response factor–dependent transcription, indicating a role of ABLIM1 in sarcomeric stress signaling (Barrientos et al., 2007). It should be noted that the potential interference with the binding of the Z1Z2 segment is just one of many mechanisms by which passive tension can be regulated. For example, there is a precedent for accessory proteins such as Ankrd1 and FHL2 to modulate passive tension by partitioning the I-band region of titin away from the thin filament (Pijl et al., 2021; Sun et al., 2024; Zhou et al., 2021). Although this partitioning mechanism appears critical for these specific proteins, it has not yet been evaluated for ABLIM1. Consequently, our hypothesis that ABLIM1 may instead function by interfering with the binding of the titin Z1Z2 segment to TCAP remains the most consistent explanation based on available data.

The IDRs within ABLIM1 are important in regulating myofilament function. One attribute is its ability to exist as an ensemble of conformations, which allows it to interact with multiple components of the myofilament, including actin and myosin, and modulate their activities in response to changing conditions, including pH and ionic strength (Sun and Kekenes-Huskey, 2023). This is because of the rugged energy landscape of ABLIM1 that enables it to sample a wide range of potential functional states. This property is particularly advantageous in the context of muscle contraction, where ABLIM1 can rapidly adapt to changing conditions to maintain optimal contractility.

PTMs provide a basis for many proteins to modulate the accessibility of binding sites and surfaces important for recognizing targets (Johnson and Lewis, 2001). Phosphorylation, in particular, alters the charge distribution of ABLIM1, shifting its conformational ensemble and influencing its interactions with other proteins, as demonstrated in our study. Supporting this are data suggesting that mutations S214A and S431A in ABLIM1 have been shown to disrupt cancer migration (Dong et al., 2020), potentially by disrupting its interactions with the actin cytoskeleton. Along these lines, we previously identified (Sun and Kekenes-Huskey, 2023) tens of PTM sites in the PhosphoSitePlus database, of which the majority were phosphorylations, but we also identified acetylation and ubiquitination sites. Similar to PTMs, mutations in IDR have the potential to impact the conformational ensemble or target binding. It is not surprising that many mutations are found in disordered regions of MAPIDs, as reported in a recent review (Sun and Kekenes-Huskey, 2023). Although missense variants have been identified in ABLIM1 in gnomAD (Chen et al., 2024) (data not shown), they remain uncharacterized and to our knowledge, a pathological variant has not been identified.

Our previous analysis of ABLIM1 (Sun and Kekenes-Huskey, 2023) revealed a large number of phosphorylation-based PTMs, many of which were in its IDRs, suggesting a complex regulatory mechanism. Coimmunoprecipitation experiments have identified numerous targets of ABLIM1, including the titin Z1Z2 segment (Stachowski-Doll et al., 2022). To expand on these hits, we queried the ELM webserver (Van Roey et al., 2014) for SLIMs in ABLIM1 and cross-referenced them with phosphorylated fragments identified in mouse data from Stachowski-Doll et al. (2022). The ELM results indicate that ABLIM1 harbors dozens of potential SLIMs, many of which are localized to its N-terminal half, of which some are at the MSSSP segment (Table S3). In particular, the MSSSP sequence is a potential target for several proteins, including GSK3β, cyclin-dependent kinase (CDK), casein kinase 2, polo-like kinase 1, and phosphoinositide-dependent protein kinase 1 (e.g., mitogen-activated protein kinases). Beyond this region, GSK3β phosphorylation motif sites are found throughout the ABLIM1 sequence. In line with these findings, knocking out GSK3β resulted in decreased, but not completely negated, phosphorylation of ABLIM1, suggesting that the altered activity of other kinases such as CDK could lead to further changes. At a minimum, these data are highly suggestive of ABLIM1 being under multiple modes of regulatory control by GSK3β and potentially other kinases.

We acknowledge key limitations and assumptions that could motivate future investigations. First, the uncertainty surrounding the interaction between LIM domains and Z1Z2 segments makes it difficult to validate our predicted docking poses. Furthermore, rigid docking may obscure the actual binding mode that would result from proteins relaxing upon binding, which could lead to inaccurate predictions. Additionally, although we suggest that ABLIM1 can bind to the Z1Z2 segment and previous studies have established a link between GSK3β KO and decreased passive tension (Stachowski-Doll et al., 2022), a direct link remains to be determined. Resolving the number of myosin heads available for recruitment using methods such as small-angle x-ray scattering may help (Irving et al., 2011) in this regard. Also, other z-repeats within titin could in principle compensate for Z1Z2 segment/TCAP interactions disrupted by ABLIM1 binding. However, the association between dilated cardiomyopathy and Z1Z2 segment mutations that reduce TCAP binding suggests that the titin Z domains do not fully compensate for Z1Z2 segment perturbation. Finally, our current atomistic simulations using phosphomimetic constructs may underestimate the structural consequences of true phosphorylation. Nonetheless, our phosphomimetic structures qualitatively reproduce phosphorylation-induced conformation changes albeit with subtle differences, which is consistent with reports for other IDRs (Yadahalli et al., 2019). Nonetheless it remains possible that more complex structural behavior like aggregation may not be recapitulated using phosphomimetics (Paleologou et al., 2008).

The key results of our study reveal significant information about the role of ABLIM1 in the sarcomere. In particular, the LIM domains of ABLIM1 were found to be capable of binding titin Z1Z2 segment, and its intrinsically disordered chain properties were shown to depend on phosphorylation. Importantly, the protein likely spans the lattice spacing in the z-disk of the sarcomere in a phosphorylation-dependent manner. These findings have broad implications for our understanding of HF. For example, the decrease in phosphorylation during HF (Jacques et al., 2008; Kirk et al., 2014) suggests that PTM control of LIM-containing proteins may play a critical role in disease progression. Together with the observed decrease in LDA in HF versus WT conditions from Stachowski-Doll et al. (2022), this highlights the importance of phosphorylation in tuning the functional contribution of ABLIM1 to cardiac contractility.

All code written in support of this publication is publicly available at https://github.com/huskeypm/pkh-lab-analyses. Simulation input files and generated data are available upon request.

Henk L. Granzier served as editor.

We are grateful for Stefano Sala’s careful reading of the manuscript.

The research reported in this publication was supported by the Maximizing Investigators’ Research Award (R35) from the National Institute of General Medical Sciences of the National Institutes of Health under grant number GM148284 and National Science Foundation (NSF) Award Numbers 2216567. This work used PSC Bridges from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support program, which is supported by NSF grant #2138259, #2138286, #2138307, #2137603, and #2138296.

Author contributions: Bin Sun: Conceptualization, data curation, methodology, software, supervision, and writing—review and editing. Alec Loftus: data curation, formal analysis, investigation, methodology, software, validation, visualization, and writing—original draft, review, and editing. Brandon Beh Goh Beh: data curation, formal analysis, investigation, methodology, visualization, and writing—original draft. Aalaythia Hepburn: investigation. Jonathan A. Kirk: conceptualization and writing—review and editing. Peter M. Kekenes-Huskey: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, and writing—original draft, review, and editing.