Identification and dynamics of the DHHC16-DHHC6 palmitoylation cascade

S-Palmitoylation is the only reversible post-translational lipid modification. Knowledge about the DHHC family of palmitoyltransferases is very limited. Here we show that mammalian DHHC6, which modifies key proteins of the endoplasmic reticulum, is controlled by an upstream palmitoyltransferase, DHHC16, revealing the first palmitoylation cascade. Combination of site specific mutagenesis of the three DHHC6 palmitoylation sites, experimental determination of kinetic parameters and data-driven mathematical modelling allowed us to obtain detailed information on the 8 differentially palmitoylated DHHC6 species. We found that species rapidly interconvert through the action of DHHC16 and the Acyl Protein Thioesterase APT2, that each species varies in terms of turnover rate and activity, altogether allowing the cell to robustly tune its DHHC6 activity.

about the DHHC family of palmitoyltransferases is very limited. Here we show that 23 mammalian DHHC6, which modifies key proteins of the endoplasmic reticulum, is 24 controlled by an upstream palmitoyltransferase, DHHC16, revealing the first 25 palmitoylation cascade. Combination of site specific mutagenesis of the three DHHC6 26 palmitoylation sites, experimental determination of kinetic parameters and data-driven 27 mathematical modelling allowed us to obtain detailed information on the 8 differentially 28 palmitoylated DHHC6 species. We found that species rapidly interconvert through the 29 action of DHHC16 and the Acyl Protein Thioesterase APT2, that each species varies in 30 terms of turnover rate and activity, altogether allowing the cell to robustly tune its 31 DHHC6 activity. auto-radiography. Decay of newly synthesized DHHC6 was biphasic, with 40% 166 undergoing gradual degradation during the first 5 hrs, and the remaining 60% 167 undergoing degradation at a greatly reduced rate (Fig. 3A). The overall apparent half-life 168 was !/! !"" ≈16hrs. To enhance DHHC6 palmitoylation, we overexpressed DHHC16 or 169 silenced APT2 expression. Both these genetic manipulations led to a dramatic 170 acceleration of DHHC6 decay, with !/! !"" ≈ 2hrs and 3hrs upon DHHC16 over expression 171 and APT2 silencing, respectively (Fig. 3A). Remarkably, mutation of Cys-328, but not of 172 Cys-329 or Cys-343, abolished the sensitivity to DHHC16 overexpression or APT2 173 silencing (Fig. 3BC, Fig. 3-figure supplement 1AB). APT2 silencing leads to 174 ubiquitination of DHHC6 (Fig. 3D) and DHHC6 degradation can be rescued by the 175 proteasome inhibitor MG132 (Fig. 3E). Thus altogether these observations indicate that 176 For DHHC6, only the non-palmitoylated form was detected, consistent with the 297 prediction (Fig. 5FG). Indeed given the dynamic range and sensitivity of Western blots, a 298 band with a ≈7-fold lower intensity than the C 000 DHHC6 band would not be detectable. 299 The model altogether predicts that under our experimental conditions, each DHHC6 300 molecule undergoes multiple rounds of palmitoylation-depalmitolyation during its life 301 cycle mostly on Cys-328. Only a subset of molecules acquire palmitate on the two other 302 sites. At any given time, approximately 60% of the molecules are not acylated. 303 304

Multi-site palmitoylation prevents large fluctuations in DHHC6 protein levels 305
Analysis of the activity of the DHHC6 mutants ( Fig. 4H) showed that the most active 306 variants are those with a cysteine at position 328, suggesting that palmitoylation at this 307 position influences palmitoyltransferase activity. However palmitoylation at this site also 308 targets DHHC6 to ERAD and thus leads to a decrease in the cellular content of DHHC6 309 (Fig. 3D). These effects raise the question of how cells can increase cellular DHHC6 310

activity. 311
We first predicted the consequences of DHHC6 hyperpalmitoylation. Calculations were 312 performed under conditions of DHHC16 overexpression and of APT2 silencing, 313 simulated by an absence of APT2. DHHC16 overexpression led to a major shift in 314 species distribution, C 011 then representing 60% of the population (Fig. 6A). Interestingly, 315 the overall DHHC6 content dropped by only 15% (Table S5). APT2 silencing in contrast 316 led to a 72% drop in DHHC6 content (Table S5). Since removing all depalmitoylating 317 activity is likely to freeze the system, we also tested a decrease of APT2 to 10% of 318 normal levels. The DHHC6 content rose by 20% and as upon DHHC16 overexpression 319 C 011 became the most abundant form (Fig. 6A). We made similar calculations for the 15 CAA mutant. Under control conditions, the total CAA content is predicted to be 84% of 321 WT, and mostly in the non-palmitoylated state (Fig. 6B). The CAA cellular content is 322 however predicted to drop to 32% of WT control levels upon DHHC16 overexpression 323 (Fig. 6B). All together, these predictions suggest that the presence of 3 palmitoylation 324 sites, as opposed to just Cys-328, renders the DHHC6 protein content more robust to 325 changes in DHHC16 activity. 326 Stochastic simulations predicts that overexpression of DHHC16 drastically increased the 327 dynamics of the network (Fig. 6C vs. 5E). Indeed, all DHHC6 molecules explored all the 328 palmitoylation states, with extremely short residence times in each (Fig. 6C). 329 Interestingly, flux analysis indicated that the abundance of C 011 was primarily due to two 330 4-step paths: both started with C 000 to C 100 , followed by palmitoylation on either of the 331 two remaining sites, subsequent depalmitoylation of Cys-328, and finally palmitoylation 332 of the remaining site ( Fig. 5-figure supplement 3). When analysing 10'000 simulations, 333 22'000 events of palmitoylation-depalmitoylation emanated from C 011 whereas only 334 4'000 events occurred between C 000 and C 100 . Thus when activity of DHHC16 is high, 335 C 011 appears to be the hub of the system. C 011 is also which is the most stable state for 336 DHHC6, explaining that the protein levels are stable. 337 We next sought to validate these predictions experimentally. Western blot analysis of 338 protein abundance showed that WT DHHC6 protein levels drastically dropped upon 339 siRNA of APT2 (Fig. 4B), but not upon DHHC16 overexpression (Fig 6DE). DHHC16 340 overexpression however led to a 60% drop in CAA expression. As expected, expression 341 of AAA was not affected (Fig. 6DE). 342 343

APEGS assay 478
The level of protein S-palmitoylation was assessed as described (Yokoi et al, 2016), with 479 minor modifications. Hela cells were lysed with the following buffer (4% SDS, 5 mM 480 EDTA, in PBS with complete inhibitor (Roche)). After centrifugation at 100,000 × g for 15 481 min, supernatant proteins were reduced with 25 mM TCEP for 1 h at 55°C or at room 482 temperature (RT), and free cysteine residues were alkylated with 20 mM NEM for 3 h at 483 RT to be blocked. After chloroform/methanol precipitation, resuspended proteins in PBS 484 with 4% SDS and 5 mM EDTA were incubated in buffer (1% SDS, 5 mM EDTA, 1 M 485 NH 2 OH, pH 7.0) for 1 h at 37°C to cleave palmitoylation thioester bonds. As a negative 486 control, 1 M Tris-HCl, pH7.0, was used. After precipitation, resuspended proteins in PBS 487 with 4% SDS were PEGylated with 20 mM mPEGs for 1 h at RT to label newly exposed 488 cysteinyl thiols. As a negative control, 20 mM NEM was used instead of mPEG (-PEG). 489 After precipitation, proteins were resuspended with SDS-sample buffer and boiled at 490 95°C for 5 min. Protein concentration was measured by BCA protein assay. HeLa cells were seeded on 12mm glass coverslips (Marienfeld GmbH, Germany) 24hr 500 prior to transfection. PAT6-myc plasmids were transfected using Fugene 501 (Promega,USA) for 48hrs. Cells were then fixed using 3% paraformaldehyde for 20min 502 at 37°C, quenched 10min with 50mM NH 4 Cl at RT and permeabilized with 0.1% TX100 503 for 5min at RT and finally blocked overnight with 0.5% BSA in PBS. Cells were washed 504 3x with PBS in between all the steps. Cells were then stained with anti-myc and anti-BiP 505 antibodies for 30min at RT, washed and incubated again 30min with their corresponding 506 fluorescent secondary antibodies. Finally cells were mounted on glass slide using 507 mowiol. Imaging was performed using a confocal microscope (LSM710, Zeiss, 508 Germany) with a 63x oil immersion objective (NA 1.4). 509  In the Goldbeter study, the model was mathematically described using mass action 518 terms. Here, due to the presence of multiple modification events, which require the 519 definition of a consistent number of parameters, we described model reactions using so-520 called "total quasi-steady state approximation" (tQSSA) (Pedersen et al, 2008). With 521 respect to the standard quasi-steady state approximation (QSSA), tQSSA is valid also 522 when the enzyme-substrate concentrations are comparable (Borghans et al, 1996).  The model can be divided in two parts; we first have a phase of synthesis of DHHC6, 529 which initially is unfolded (U in Fig. 5A) and subsequently reaches the folded initially 530 non-palmitoylated (C 000 ) form. Each site can subsequently be palmitoylated leading to 531 C 100 C 010 and C 001 . These species can then undergo a second palmitoylation step 532 leading to C 110 , C 101 and C 011 , and a third leading to C 111 . Since palmitoylation is 533 reversible, palmitate can be removed from each of the sites, from all of the species. • The three palmitoylation sites may have different affinities with respect to the 544 palmitoylation/depalmitoylation enzymes, therefore we adopted separated Kms 545 and V max s for the different sites. 546 • All the palmitoylation steps are reversible. APT catalyses the depalmitoylation 547 steps. 548 • Palmitate was considered to be available in excess. 549 A detailed account of all reactions and differential equations in the model is given in 550 Tables S1-3. The types of data used, the description of the Genetic Algorithm (GA) and the step-by-559 step application of the optimization algorithm used to find values for the parameters are 560 described in detail in (Dallavilla et al, 2016). In this paper we use the exact same 561 procedure. 562 For the parameter estimation we used a calibration set that correspond to the data 563 visible in (Fig. 5B and Fig.S5). The remaining part of the data ( Fig.5B and S6) was used 564 to validate the output of the model and to verify the accuracy of its prediction capabilities. 565 Because of the presence of multiple objectives there does not exist a single solution that 566 simultaneously optimizes each objective, so the algorithm provides as output a local 567 Pareto set of solutions, which are equally optimal with respect to the fitness function we 568 defined. From the Pareto set provided by the GA and in order to be more accurate and 569 to reduce the variability in the output of the model, we selected the set of parameters 570 that best fitted the calibration data. The procedure for the selection of the subset is 571 The results of the optimization can be found in Table S3. has no detectable activity. Only upon modification by an upstream palmitoylatransferase, 606 which we identified as DHHC16, does it acquire significant transferase activity. DHHC16, 607 also known as Ablphilin 2 (Aph2), is, as DHHC6, expressed in many human tissues in 608 particular in the heart, pancreas, liver, skeletal muscle (Zhang et al, 2006). Its name 609 Ablphilin stems from its ability to interact with the non-receptor tyrosine kinase c-Abl at Our findings raise the possibility that some of these effects may involve DHHC6. 615 We found that DHHC6 palmitoylation is highly dynamic, involving Acyl Protein 616 Thioesterase 2. APT2, which was recently found to be involved in cell polarity-mediated 617  PEGylation analysis of mouse tissues indicate that in vivo DHHC6 palmitoylation is far 643 more pronounced than in tissue cultured cells (Fig. 6H). This situation was mimicked 644 experimentally by overexpression of DHHC16 and revealed the importance of having 3 645 palmitoylation sites, rather than just Cys-328. We indeed found that in the presence of 3 646 sites, the cellular DHHC6 activity could be enhanced by DHHC16 overexpression. 647 Increase DHHC16 activity led to a shift of species distribution C 011 , the most stable form, 648 becoming the most abundant. Altogether this study indicates that DHHC6 can adopt in 3 649 types of states: 1) stable and inactive (C 000 ); 2) highly active, short-lived and rapidly 650 turned over (C 100 ); 3) moderately active, long lived and highly stable (C 010 and C 011 ). 651 Such a regulatory system allows the cell to tightly control the activity of DHHC6. Why 652 excessive DHHC6 activity is detrimental to a cell or an organism will require further 653

studies. 654
By combining experimentation and the predictive power of data-driven mathematical 655 modelling, we were able to obtain unprecedented insight into the dynamics of 656 palmitoylation and the role and properties of single acylated species. We indeed found 657 that depending on the site of modification both the activity and turnover of DHHC6 were 658 strongly affected, indicating that palmitoylation has a more subtle and precise effect that 659 merely increasing hydrophobicity, consistent with recent findings on single acylated    DHHC6 is not autopalmitoylated. Hela cells silenced with control lentiviruses or with shDHHC6 lentiviruses were transfected with plasmids encoding WT myc-DHHC6 or mutant myc-DHHC6. Cells were then metabolically labeled 2 hours at 37°C with 3H-palmitic acid. Proteins were extracted and immunoprecipitated with anti-myc antibodies, subjected to SDS-PAGE and analyzed by autoradiography ( 3 H-palm), quantified using the Typhoon Imager or by immunoblotting with anti-myc antibodies. 3H-palmitic acid incorporation into DHHC6 values were calculated and were set to 100% for WT DHHC6 constructs and all DHHC6 mutants were expressed relative to this. N=4. B. All DHHC6 cysteine are palmitoylated by DHHC16. Hela cells were transfected with siRNA silencing indicated DHHC for 72h and with WT and cysteine mutants myc-tagged DHHC6 construct for the last 24h. Cells were then metabolically labeled 2 hours at 37°C with 3H-palmitic acid. Proteins were extracted, immunoprecipitated with myc antibodies and subjected to SDS-PAGE and analyzed by autoradiography, quantified using the Typhoon Imager or by immunoblotting with myc antibodies. 3 H-palmitic acid incorporation into WT and mutants DHHC6 constructs were quantified and normalized to protein expression level. The calculated value of 3 H-palmitic acid incorporation into WT or mutants DHHC6 was set to 100% for a non relevant siRNA (Ctrl) and all siRNA were expressed relative to this. N=3. C. Palmitoylation of WT DHHC6 after human DHHC overexpression. Hela cells were transfected with indicated DHHC constructs for 24h with WT myc-tagged DHHC6 construct. Cells were then metabolically labeled 2 hours at 37°C with 3H-palmitic acid. Proteins were extracted, immunoprecipitated with myc antibodies and subjected to SDS-PAGE and analyzed by autoradiography, quantified using the Typhoon Imager or by immunoblotting with myc antibodies. 3 H-palmitic acid incorporation into WT DHHC6 constructs were quantified and normalized to protein expression level. The calculated value of 3H-palmitic acid incorporation into WT DHHC6 was set to 100% for a non relevant plasmid (Ctrl) and all DHHC were expressed relative to this. N=5. D. Co-immunoprecipitation of DHHC6 with DHHC16. Hela cells were transfected with plasmids encoding WT myc-DHHC6 and flag-tagged DHHC16 constructs for 24h. Proteins were extracted, a total cell extract was analyzed (TCE) and proteins were immunoprecipitated with myc or flag antibodies and subjected to SDS-PAGE, then analyzed by immunoblotting with anti-myc or anti-flag antibodies.  B.DHHC6 is palmitoylated but not DHHC16. Hela cells were transfected with WT myc-tagged DHHC6 or myc-tagged DHHC16 construct for 24h. Cells were then metabolically labeled 2 hours at 37°C with 3H-palmitic acid. Proteins were extracted, immunoprecipitated with myc antibodies and subjected to SDS-PAGE and analyzed by autoradiography, quantified using the Typhoon Imager or by immunoblotting with myc or flag antibodies. C. Analysis of protein acylation in Hela cells (Ctrl). Hela were transfected 24h with WT myc-DHHC6 or WT flag-DHHC16 constructs. Cell membranes were recovered by centrifugation and incubated with MMTS and then with hydroxylamine (+HA) or with TRIS (-HA) together with free thiol group binding beads. Eluted fractions were analyzed by immunoblotting with the indicated antibodies. The input fraction (PNS) was loaded as 1/10 this amount. DHHC6 were immunoprecipitated and subjected to SDS-PAGE and analyzed by autoradiography, quantified using the Typhoon Imager, and western blotting with anti-myc antibodies. 35 S-methionin/cysteine incorporation into different DHHC6 constructs were quantified for each times, normalized to protein expression level. The calculated value of 35 S-methionin/cysteine incorporation into DHHC6 was set to 100% for t=0 after the 20min pulse and all different times of chase with complete medium were expressed relative to this. N=3.
3. Degradation Where the competition terms f and b in S1 Table are defined as: Km b3 where k catP16c1 is the catalytic constant for DHHC16 when catalyzing palmitoylation of C 100 on DHHC6, while k catAPTc1 is the catalytic constant for APT2 for the same site. k c1P16b and k c1P16ub are the binding and unbinding rate constants of DHHC6 to DHHC16 on site C 100 respectively. k c1APTb and k c1APTub k c1APTub are the binding and unbinding rate constants of DHHC6 to APT1 on site C 100 respectively. In our model we assume that the binding and unbinding rate to DHHC6 and APT can be different for each state of the model.  Table S2. Mass balance equations. The following table describes the mas balance for each of the species of DHHC6 model. The rates of the mass bal each state are described in detail in S1 Table. Mass balance Mass balance Table S3. Model parameters. The output of GA is a set of optimal solutions, where a solution is a complete set of parameter needed to perform model simulations. From this set we extracted a sub-set of 152 solutions that obtained a GA score better than a set threshold for each objective. During the analysis the model was simulated for each set of parameters of the sub-set. We then reported in this paper the mean of the outputs along with the 1 st and 3 rd quartile of their distribution.