Raphaëlle Versini

Marc Baaden

Laetitia Cavellini

Mickaël M Cohen

Antoine Taly
email: taly@ipbc.fr

Patrick F J Fuchs
email: patrick.fuchs@u-paris.fr

Lys716 in the transmembrane domain of yeast mitofusin Fzo1 modulates anchoring and fusion

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Mitochondria form a complex network inside the cells, undergoing continuous fusion and fission events. These processes shape mitochondrial dynamics and are essential for the maintenance, function, distribution and inheritance of mitochondria. The morphology of the latter therefore respond to the everchanging physiological changes of the cell [START_REF] Westermann | Mitochondrial fusion and fission in cell life and death[END_REF].

Mitofusins are large GTPase transmembrane proteins involved in the tethering to and fusion of mitochondrial outer membranes (OM) [START_REF] Ramachandran | Mitochondrial dynamics: The dynamin superfamily and execution by collusion[END_REF]. The mitofusins Mfn1 and Mfn2 are found in mammals [START_REF] Ishihara | Mitofusin 1 and 2 play distinct roles in mitochondrial fusion reactions via GTPase activity[END_REF][START_REF] Rojo | Membrane topology and mitochondrial targeting of mitofusins, ubiquitous mammalian homologs of the transmembrane GTPase fzo[END_REF]. Fzo1 (Fuzzy Onion 1 ) is the only mitofusin homolog in Saccharomyces cerevisiae [START_REF] Hermann | Mitochondrial fusion in yeast requires the transmembrane GTPase fzo1p[END_REF]. The structure of Mfn1 was partially solved, but without its transmembrane domain [START_REF] Cao | Mfn1 structures reveal nucleotidetriggered dimerization critical for mitochondrial fusion[END_REF][START_REF] Yan | Structural basis for gtp hydrolysis and conformational change of mfn1 in mediated membrane fusion[END_REF], and no solved structures are available for either Mfn2 or Fzo1. Mitochondrial inner membrane fusion and cristea organisation are mediated by human OPA1 (Optic Atrophy 1 ) [START_REF] Wong | The dynamin-related GTPase, mgm1p, is an intermembrane space protein required for maintenance of fusion competent mitochondria[END_REF] and yeast Mgm1 (Mitochondrial Genome Maintenance 1 ) [START_REF] Frezza | OPA1 controls apoptotic cristae remodeling independently from mitochondrial fusion[END_REF].

Mitochondrial fusion dysfunction can cause neurodegenerative diseases such as Parkinson's, Alzheimer's, and Huntington's [START_REF] Chen | Mitochondrial dynamics-fusion, fission, movement, and mitophagy-in neurodegenerative diseases[END_REF][START_REF] Knott | Mitochondrial fragmentation in neurodegeneration[END_REF] as well as cancer [START_REF] Keshav | Special issue: Mitochondria in cancer[END_REF][START_REF] Allegra | Relationship between mitofusin 2 and cancer[END_REF]. A number of studies have linked Mfn1 and Mfn2 to various cancers, including breast cancer [START_REF] Darvish Moghaddam | Apoptotic effects of melittin on 4t1 breast cancer cell line is associated with up regulation of mfn1 and drp1 mrna expression[END_REF], liver cancer [START_REF] Zhang | Mfn1-dependent alteration of mitochondrial dynamics drives hepatocellular carcinoma metastasis by glucose metabolic reprogramming[END_REF], lung cancer [START_REF] Jalees | Inhibition of mitochondrial fission prevents cell cycle progression in lung cancer[END_REF], cervical cancer [START_REF] Ahn | Mitofusin-2 promotes the epithelial-mesenchymal transition-induced cervical cancer progression[END_REF], and colon cancer [START_REF] Cheng | Prognostic impact of mitofusin 2 expression in colon cancer[END_REF]. They showed that their dysregulated expression was associated with increased cell proliferation, invasion, and resistance to chemotherapy. Regulation of mitofusin activity has been shown in preclinical studies to reduce cancer cell growth and spread. In addition, research has shown that mutations in Mfn2 trigger the development and progression of muscular dystrophies such as Charcot-Marie-Tooth type 2A, the most common form of axonal CMT disease [START_REF] Feely | Mfn2 mutations cause severe phenotypes in most patients with cmt2a[END_REF][START_REF] Stephan | Mutations in the mitochondrial gtpase mitofusin 2 cause charcot-marie-tooth neuropathy type 2a[END_REF]. The precise mechanism by which mitofusins contribute to cancer or neurodegenerative disease, as well as the exact molecular fusion mechanism mediated by mitofusins, require further structural investigation.

From other membrane fusion machines such as the SNARE proteins, we know that anchoring in the membrane, transmembrane domain amino-acid composition, and dynamic properties are crucial elements that determine the interplay between the main players of the fusion process at different stages [START_REF] Lindau | Coarse-Grain Simulations Reveal Movement of the Synaptobrevin C-Terminus in Response to Piconewton Forces[END_REF][START_REF] Han | Synaptobrevin transmembrane domain determines the structure and dynamics of the snare motif and the linker region[END_REF][START_REF] Wehland | Role of the transmembrane domain in snare protein mediated membrane fusion: peptide nucleic acid/peptide model systems[END_REF][START_REF] Dhara | v-snare transmembrane domains function as catalysts for vesicle fusion[END_REF]. With the help of a membrane-embedded structural model of Fzo1, such properties, the fusion itself, as well as its intermediate stages, could be studied by computational approaches, in particular molecular dynamics simulations. To this end, we are striving to build a complete structural model of Fzo1, including in particular the fusion-related transmembrane domains, for which we have very limited data.

Fzo1 has two heptad repeat domains (HRs) with coiled-coils on its N-terminal side: HRN (only in yeast) and HR1 flanking a GTPase domain. A third heptad repeat domain HR2 is located at the Cterminus [START_REF] Hermann | Mitochondrial fusion in yeast requires the transmembrane gtpase fzo1p[END_REF][START_REF] Griffin | Domain interactions within fzo1 oligomers are essential for mitochondrial fusion[END_REF]. Some models of Fzo1 are based on the mitofusin-related bacterial dynamin-like protein (BDLP) [START_REF] Harry | A bacterial dynamin-like protein[END_REF][START_REF] Harry | Structure of a bacterial dynaminlike protein lipid tube provides a mechanism for assembly and membrane curving[END_REF]. BDLP is involved in membrane remodelling and exists in two conformational states, a closed, compact version that transitions to an open, extended structure upon GTP binding. The latter structure can therefore bind lipid bilayers and results in self-accelerating polymerization leading to a coated lipid bilayer and strong curvature. However, BDLP has no transmembrane part, whereas Fzo1 is embedded in the mitochondrial OM via two transmembrane segments exposing N-and C-terminal parts to the cytosol and a loop to the intermembrane space [START_REF] Fritz | Connection of the mitochondrial outer and inner membranes by fzo1 is critical for organellar fusion[END_REF]. The two transmembrane segments are putative α-helices called TM1 and TM2 [START_REF] Vecchis | A membrane-inserted structural model of the yeast mitofusin fzo1[END_REF].

We previously constructed an initial model of Fzo1 based on the closed compact conformation of BDLP, which was experimentally validated by mutation studies [START_REF] Vecchis | A membrane-inserted structural model of the yeast mitofusin fzo1[END_REF]. The opened extended structure of BDLP was used to create a second model [START_REF] Brandner | Physicsbased oligomeric models of the yeast mitofusin fzo1 at the molecular scale in the context of membrane docking[END_REF][START_REF] Vecchis | A molecular perspective on mitochondrial membrane fusion: From the key players to oligomerization and tethering of mitofusin[END_REF]. In the absence of a BDLP transmembrane template, the structure of TM1 and TM2 of Fzo1 was predicted using an ab initio method, namely the PREDDIMER [START_REF] Polyansky | Preddimer: a web server for prediction of transmemrane helical dimers[END_REF] web server. Here we attempt a more comprehensive investigation of the possible membrane associations of the two helices to form a dimer using coarse-grained molecular dynamics simulations. Indeed, a new version of the popular Martini force field solved a previous bias favoring protein/protein association in membranes [START_REF] Paulo | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF]. This new force field, Martini 3, therefore paved the way to more realistic studies of membrane protein associations [START_REF] Amita R Sahoo | Transmembrane dimers of type 1 receptors sample alternate configurations: Md simulations using coarse grain martini 3 versus alphafold2 multimer[END_REF]. In the case of Fzo1, the TM domain presents a particular challenge, as there is a basic conserved lysine residue in the middle of TM1. In addition to dimerization itself, we are investigating the extent to which the protonation state of this residue may affect dimer stability, membrane incorporation, and destabilization. Lysine protonation may be functionally relevant, as the literature has hypothesized specific roles for individual transmembrane domains in fusion and, in particular, charge-related effects. Specifically for the SNARE system, Lindau et al. postulated a mechanism in which the formation of fusion pores is triggered by the movement of the charged C-terminus of the transmembrane helix syb2 within the membrane [START_REF] Lindau | Coarse-Grain Simulations Reveal Movement of the Synaptobrevin C-Terminus in Response to Piconewton Forces[END_REF].

Through extensive molecular simulations across coarse-grained and atomistic scales, we arrive at an improved model of the Fzo1 transmembrane domain that reveals new insights into the dimerization interface, key regulatory residues, and their potential influences on mitochondrial fusion.

Materials and Methods

The transmembrane (TM) domain of Fzo1 contains two TM α-helices (called TM1 and TM2) connected by a loop, plus some flanking residues on the N-terminal part of TM1 and C-terminal part of TM2. In 2/31 addition, TM1 possesses a lysine residue (Lys716) at its exact center, thus located in the middle of the membrane. This location therefore raises the question of its protonation state, knowing that the local pKa of lysine may go below 7.0 in very apolar environments [START_REF] Maccallum | Distribution of amino acids in a lipid bilayer from computer simulations[END_REF][START_REF] Gleason | Buried lysine, but not arginine, titrates and alters transmembrane helix tilt[END_REF][START_REF] Panahi | Membrane environment modulates the pka values of transmembrane helices[END_REF].

The computational protocol used in this study consisted in two main parts. The first part used coarsegrained (CG) simulations for sampling the TM1/TM2 contacts and selecting the best possible TM1/TM2 dimer. Two CG force-fields (Martini 2 and Martini 3) were tested as well as the two protonation states of Lys716 (neutral or positively charged). In the second part, the best CG dimer was chosen, the loop between TM1 and TM2 added, the structure backmapped to an all-atom (AA) representation, followed by AA simulations to refine the whole Fzo1 TM domain. Everything is summarized in Figure 1. Below, we explain the details of each part. The GROMACS 2018.5 program was used to perform all molecular dynamics (MD) simulations [START_REF] James | Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF].

Coarse-Grained conformational sampling

Peptide Definition

The S. cerevisiae Fzo1 sequence and domains were obtained from the UniProtKB database (Universal Protein Resource Knowledgebase [START_REF]UniProt: a worldwide hub of protein knowledge[END_REF], entry: P38297). The secondary structures of the protein were predicted using PSIPRED 4.0 [START_REF] Daniel | The PSIPRED protein analysis workbench: 20 years on[END_REF][START_REF] David | Protein secondary structure prediction based on position-specific scoring matrices[END_REF]. The positions of the transmembrane helices of Fzo1 were determined using TMHMM2.0 [START_REF] Krogh | Predicting transmembrane protein topology with a hidden Markov model: application to complete genomes[END_REF][START_REF] Sonnhammer | A hidden markov model for predicting transmembrane helices in protein sequences[END_REF]. Some peptides were defined (PEPT1 or PEPT2) which contained either TM1 or TM2 respectively, plus various flanking residues: one from Arg695 to Lys735 called PEPT1-41 (of 41 residues length), another from Ser702 to Ser730 called PEPT1-29 (of 29 residues length), and a last one from Ser733 to Leu761 called PEPT2 (of 29 residues length). PEPT1-29 and 3/31 PEPT2 were used with Martini 2, whereas PEPT1-41 and PEPT2 were used with Martini 3. All sequences are shown in Figure 1A. The peptide 3D-structures were built using Basic Builder 1.0.2 on the Mobyle server (https://mobyle.rpbs.univ-paris-diderot.fr) [START_REF] Néron | Mobyle: a new full web bioinformatics framework[END_REF]. The residues from Thr706 to Ile726 (TM1) and from Leu737 to Leu57 (TM2) were defined as a perfect α-helix and the flanking residues were defined as coil. This step gave a first all-atom 3D structure for each peptide, which was then used as an input to Martini and CHARMM-GUI tools described in the next section.

Box setup

When we started the work, Martini 2 was the main CG Martini force field and Martini 3 was in its beta version. However, the Martini 3 force field came out during the course of this work [START_REF] Paulo | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF]. We thus decided to test both versions of the force-field, Martini 2 and Martini 3. We tested the two possible charge-states of Lys716 (positively charged or neutral), resulting in a total of four systems.

The Martini 2 systems were built using the CHARMM-GUI Membrane builder [START_REF] Jo | Charmm-gui: A web-based graphical user interface for charmm[END_REF][START_REF] Qi | Charmmgui martini maker for coarse-grained simulations with the martini force field[END_REF][START_REF] Djurre | Improved parameters for the martini coarse-grained protein force field[END_REF] with the Martini 2 coarse-grained (CG) force field [START_REF] Siewert J Marrink | The martini force field: coarse grained model for biomolecular simulations[END_REF][START_REF] Monticelli | The martini coarse-grained force field: extension to proteins[END_REF]. The two peptides were separately embedded in a mixed lipid bilayer of 50 palmitoyl-oleoyl-phosphatidylcholine (POPC) and 50 palmitoyl-oleoylphosphatidylethanolamine (POPE) (these two lipids are the most abundant in mitochodrion OM). Thus we had one box with PEPT1 (PEPT1-29 or PEPT1-41) and another one with PEPT2. About 40 water molecules per lipid were used as solvent for both PEPT1-29 and PEPT2, and about 50 water molecules per lipids were placed in the system of PEPT1-41. To the solvent were added 0.15M of NaCl for each system. The two systems were then assembled, in order to have the two peptides embedded in the same membrane at a distance of about 5.2 nm. In order to obtain a neutral Lys716 system, we replaced in the ITP file the Qd charged bead by a P1 neutral and polar bead.

Since the mapping between Martini 3 and Martini 2 is different, we could not convert the topologies directly. Instead, the Martini 3 systems were built from atomistic models of the peptides (using the same coordinates than for Martini 2) using the martinize2.py workflow module (see https://github. com/marrink-lab/vermouth-martinize). These Martini 3 structures of the peptides were then placed in lieu of the peptide structures in the Martini 2 system. In order to obtain a neutral Lys716 system, we replaced in the ITP file the SQ4p charged bead by a SN6d neutral and polar bead. As the mapping of the lipids did not change between the two versions of the force-field, we were able to use the membranes produced with Martini 2.

Simulation parameters

For all systems, we followed the protocol of CHARMM-GUI which consisted in two energy minimizations of 5000 steps, followed by an equilibration of 5 simulations within the NPT ensemble. The velocityrescaling thermostat [START_REF] Bussia | Canonical sampling through velocity rescaling[END_REF] at 303.15 K and the Berendsen barostat at 1 bar [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] were applied for a sequence of 5 simulations of 1 ns, 1 ns, 1 ns, 0.75 ns and 1 ns (with a timestep of 0.002 ps, 0.005 ps, 0.01 ps, 0.015 ps, 0.020 ps respectively). From the equilibrated system, a last step of 1 ns equilibration (time step 0.020 ps) with the same thermostat and barostat parameters was then used to redefine the velocities of the system (with a fixed seed). A production run of 10 µs followed, at 303.15 K using the velocity-rescaling thermostat [START_REF] Bussia | Canonical sampling through velocity rescaling[END_REF] (lipid, water and proteins coupled separately) and 1 bar using the Parrinello-Rahman barostat [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] (compressibility of 3.0 × 10 -4 bar -1 ). Pressure coupling was applied semi-isotropically. These final two steps were repeated multiple times, the seed being changed for each preceding equilibration step. A time step of 0.020 ps was used with the leapfrog integrator. Lennard-Jones interactions were cutoff at 1.1 nm. Bond lengths were constrained using the LINCS algorithm [START_REF] Hess | LINCS: A linear constraint solver for molecular simulations[END_REF]. The reaction-field method [START_REF] Ilario | A generalized reaction field method for molecular dynamics simulations[END_REF] was used for evaluating electrostatic interactions, with a Coulomb distance cutoff of 1.1 nm, a relative dielectric constant of 15. The neighbor list was updated every 20 steps.

4/31

Clustering and model extraction

We ran as many simulations as necessary in order two obtain 25 trajectories ending in successful dimerization for each of the four conditions (two Martini versions and two protonation states). The unsuccessful trajectories were discarded. The successful ones were submitted to a conformational based clustering. The first 2 µs of the productions were systematically ignored. All simulations were then concatenated (using 161 frames per simulation representing 8 µs) resulting in a total of 4025 frames. This concatenated trajectory was used as input to the gmx cluster program of GROMACS. The GROMOS method [START_REF] Daura | Peptide folding: When simulation meets experiment[END_REF] was used. Briefly, one first computes the pairwise RMSD matrix of TM1 and TM2 (RMSD of all pairs of conformations, considering backbone beads only) and then group conformations in different clusters. Within each cluster, any pair of conformations presents an RMSD below a cutoff that has to be chosen. This non-supervised method (the number of clusters is not fixed and is a result of the clustering) maximizes the size of the clusters. The cutoff was set to 0.3 nm in order to get 5 to 6 clusters totalizing 80 to 90 % of all conformations.

Trajectory analysis

The crossing angle between the TM domains was calculated as described by Chothia et al. [START_REF] Chothia | Helix to helix packing in proteins[END_REF] using an in-house script. Only the backbone beads of TM1 and TM2 were considered.

Contacts between TM1 and TM2 were studied using the python library MDAnalysis [START_REF] Michaud-Agrawal | MD-Analysis: A toolkit for the analysis of molecular dynamics simulations[END_REF][START_REF] Gowers | MDAnalysis: A python package for the rapid analysis of molecular dynamics simulations[END_REF]. The minimum distance between each residue pair was then plotted using the R library ggplot2 [START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF].

A principal component analysis (PCA) of the cartesian coordinates was carried out using the programs gmx covar and gmx anaeig of the GROMACS package. It was done on a combined trajectory of the dimerized structures of the 4 systems using the backbone beads of TM1 and TM2. The results are presented in terms of free energy projection against the two first principal components (PC1 and PC2) using a grid: ∆G(i, j) = -RT ln p ij where i and j are indices of the grid points, p ij the probability of finding a conformation in grid point (i, j) and ∆G(i, j) the relative free energy of that grid point. The lowest free energy was set to 0 kJmol -1 . 150 grid points were used for PC1 and PC2.

Lipid packing defects are small apolar areas of the membrane which are accessible to water. They were quantified using PackMem [START_REF] Vamparys | Conical lipids in flat bilayers induce packing defects similar to that induced by positive curvature[END_REF][START_REF] Gautier | Packmem: A versatile tool to compute and visualize interfacial packing defects in lipid bilayers[END_REF]. A protrusion event is defined as one of the carbon atoms of a lipid tail bulging into the polar layer (or a backbone bead of a lipid tail), extending 0.1 nm above (or below, depending on the leaflet) its phosphorus atom [START_REF] Larsson | Lipid tail protrusion in simulations predicts fusogenic activity of influenza fusion peptide mutants and conformational models[END_REF][START_REF] Mukarram | Solvent-exposed lipid tail protrusions depend on lipid membrane composition and curvature[END_REF]. Protrusions were identified using an in-house script. Both analyses were performed separatly for each leaflet.

All-Atom refinement of the model

The CHARMM36m force field for proteins [START_REF] Huang | CHARMM36m: An improved force field for folded and intrinsically disordered proteins[END_REF] and CHARMM36 for lipids [START_REF] Klauda | Update of the charmm all-atom additive force field for lipids: Validation on six lipid types[END_REF] were used for the remaining all-atom simulation described in this section.

System building and equilibration

Prior to all-atom simulations, the loop between TM1 and TM2 was reconstructed and sampled using CG simulations. We used position-restraints on TM1 and TM2 in order to maintain the contacts predicted in the previous phase. The details of this step can be found in the Supplementary Material. Briefly, the best loop conformation was chosen based on conformational clustering.

The whole TM domain including the loop was then backmapped using the CHARMM-GUI backmapping tool [START_REF] Tsjerk | Going backward: A flexible geometric approach to reverse transformation from coarse grained to atomistic models[END_REF]. We reduced the box size in order to have 100 lipids located around the TM domains, as well as 40 water molecules per lipids and 0.15M of NaCl. In total, the system consisted in 26144 atoms.

The preceding step introduced a layer of vacuum on the box edges. Using 8 minimizations of 5000 steps, we progressively shrinked the box size to get rid of this vacuum layer and recover periodic boundary 5/31 conditions. The system was then submitted to a sequence of 6 equilibrations of 125 ps, 125 ps, 2 ns, 2 ns, 2 ns, and 2 ns (with a timestep of 0.001 ps, 0.001 ps, 0.001 ps, 0.002 ps, 0.002 ps, 0.002 ps respectively) in which we progressively released the position-restraints on the protein. The first two equilibrations were performed within the NVT ensemble, with the Berendsen thermostat, and the following simulations were performed within the NPT ensemble, with the Berendsen thermostat and barostat. The temperature was maintained at 303.15 K and the pressure at 1 bar. Pressure was applied semi-isotropically. Electrostatic interactions were calculated with the particle-mesh-Ewald (PME) method [START_REF] Darden | Particle mesh ewald: An n log (n) method for ewald sums in large systems[END_REF][START_REF] Essmann | A smooth particle mesh ewald method[END_REF], with a real-space cutoff of 1 nm. Van der Waals interactions were computed using a Lennard-Jones force-switching function over 10 to 12 Å. Bond lengths were constrained using the LINCS algorithm [START_REF] Hess | Lincs: A linear constraint solver for molecular simulations[END_REF]. Water molecules were kept rigid with the SETTLE algorithm [START_REF] Miyamoto | Settle: An analytical version of the shake and rattle algorithm for rigid water models[END_REF].

Temperature replica-exchange molecular dynamics

In order to explore efficiently the conformational landscape of the TM domain, we used temperature replica-exchange molecular dynamics (T-REMD) [72]. The replica temperatures were predicted using the webserver https://virtualchemistry.org/remd-temperature-generator/ [START_REF] Patriksson | A temperature predictor for parallel tempering simulations[END_REF] by setting the exchange probability to 0.2 and the temperature range between 303.15 and 399.83 K. The upper limit of nearly 400 K was chosen so that conformational sampling of the protein was more efficient, but the membrane stayed intact (if heated too much, it can explode). In total, we obtained 38 replica.

Each replica was equilibrated at the chosen temperature for 1 ns (time step 0.002 ps) using the Berendsen thermostat and barostat (at 1 bar) with different starting velocities. The production run of 500 ns followed with the velocity-rescaling thermostat [START_REF] Bussia | Canonical sampling through velocity rescaling[END_REF] (lipid, water and proteins coupled separately) and the Parrinello-Rahman barostat at 1 bar [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] (compressibility of 4.5 × 10 -5 bar -1 ). Exchanges between neighboring replicas were attempted every 10 ps. The other settings were identical to those described in the previous section.

Trajectory analysis

All the analyses were performed either based on individual replica that diffuse in temperature, or on the conformations of the bottom temperature (303.15 K). On the latter, crossing-angles (using helical Cα atoms only), conformational clustering (using a cutoff of 0.2 nm and a distance between the TM1 and TM2 center of mass below 2 nm) and contact-maps were calculated using the same programs as for the CG simulations. The cutoff of 0.2 nm used for the clustering was chosen so that each cluster roughly coincided with the group of points obtained in the PCA analysis below.

A PCA was carried out on a trajectory combining the REMD structures used for the clustering and the Martini 3 neutral Lys716 structures. To do this, the all-atom conformations were converted to a CG representation, and concatenated to the Martini 3 structures for the covariance matrix calculation. Only the backbone beads of TM1 and TM2 were considered.

A conformational clustering was performed on the Cα atoms of TM1 and TM2 using the same algorithm and tool as for CG simulations (cutoff of 0.2 nm). The central structure of the first and third clusters were then extracted. We selected these two structures because of their positioning within the deepest energy wells identified by the previous PCA.

Hydrogen bonds were calculated using gmx hbond from GROMACS 2018.5. Contacts within a distance cut-off of 0.35 nm and up to 30 degree off-axis angle were considered.

AlphaFold2 predictions

To predict the structures of Fzo1, we used Alphafold version 2.2 [START_REF] Jumper | Highly accurate protein structure prediction with alphafold[END_REF][START_REF] Evans | Protein complex prediction with alphafold-multimer[END_REF], and Colabfold 1.3.0 [START_REF] Mirdita | Colabfold: Making protein folding accessible to all[END_REF], both monomer and multimer versions. The provided sequences of Fzo1 (from residues 491 to 813, total length of 323 residues, Uniprot: P38297) and Ugo1 (Uniprot: Q03327) were taken from Uniprot [START_REF]UniProt: a worldwide hub of protein knowledge[END_REF].

6/31

Plasmids,yeast strains and growth conditions.

Mutations on K716 (plasmids MC430, MC437, MC587 and MC589, Table1) were generated by PCR using QuikChange Lightning Site-Directed Mutagenesis Kit (Agilent Technologies 210518, Santa Clara, California, USA). The QuikChange Primer Design Program is available online at https://www.agilent.com/store/primerDesignProgram.jsp and was used to design mutagenic primers based on the Fzo1 sequence. Standard methods were used for growth, transformation and genetic manipulation of S. cerevisiae. Minimal synthetic media [Difco yeast nitrogen base 291940 (Voigt Global Distribution, Inc., Lawrence, Kansas, USA), and drop-out solution] supplemented with 2% dextrose (SD) or 2% glycerol (SG) were prepared as previously described in Sherman et al, J. Methods in Yeast Genetics (1986) [START_REF] Sherman | Methods in yeast genetics: a laboratory course manual[END_REF]. All experiments were performed with the fzo1D shuffle strain MCY571 where mitochondrial fusion efficiency is maintained by the pRS416-FZO1 shuffle plasmid carrying a copy of the Fzo1 wild-type gene and the URA3 selection marker. This strain was transformed with the pRS314 FZO1 plasmids described in Table 1 carrying the TRP1 selection marker and resulting transformants were plated on SD selective media lacking uracil and tryptophan. Colonies were isolated on SD selective media and replica-plated on 5'-fluoroorotic acid (5'-FOA) plates to cure the strains from the FZO1 shuffling plasmid.

Spot assays.

Cultures grown overnight in SD selective media lacking tryptophan were pelleted, resuspended at OD600=1 and serially diluted (1:10) four times in water. Five microliters of the dilutions were spotted onto SD and SG selective plates and grown for 2 days at 23 °C, 30 °C or 37 °C.

Protein extracts and immunoblotting.

Cells grown in SD selective media lacking tryptophan were collected during the exponential growth phase (OD600 = 0.7-1). Total protein extracts were prepared using the NaOH/trichloroacetic acid (TCA) lysis technique [START_REF] Volland | Endocytosis and degradation of the yeast uracil permease under adverse conditions[END_REF]. Proteins were separated on 8% SDS-PAGE gels and transferred to nitrocellulose membranes (Amersham Protran 0,45 µm 10600002; GE Healthcare, Little Chalfont, United Kingdom). The primary antibodies used for immunoblotting were monoclonal anti-Pgk1 (Abcam ab113687, Cambridge,United Kingdom) and polyclonal anti-Fzo1 (generated by Covalab, Bron, France). Primary antibodies were detected using horseradish peroxidase-conjugated secondary anti-mouse or anti-rabbit antibodies (HRP, Sigma-Aldrich 12-348 and A5278, Saint-Louis, Missouri, USA), followed by incubation with a Clarity Western ECL Substrate (Bio-Rad 1705060, Hercules, California, USA). Images of the immunoblots were acquired using a Gel DocTM XR+ (Bio-Rad) and analysed using the Image Lab 3.0.1 software (Bio-Rad).

Results

The association of transmembrane helices is reproducible with two coarsegrained force fields

Our aim was to predict the structure of the transmembrane domain of Fzo1 in a 1:1 POPC-POPE bilayer. This lipid composition was chosen considering the two most abundant phospholipids in the outer mitochondrial membrane and was the same as in our previous study [START_REF] Vecchis | A membrane-inserted structural model of the yeast mitofusin fzo1[END_REF]. To evaluate the possible associations between the two helices TM1 and TM2 in coarse-grained (CG) simulations, the domain was cut into two parts and some peptides were defined containing either TM1 or TM2: PEPT1-29, PEPT1-41, and PEPT2. Each sequence contains also some flanking residues on both the N-and Cterminal sides (see Figure 1A). In the simulations performed with the Martini 2 force field, PEPT1-29 and PEPT2 were used. However, in the simulations performed with the Martini 3 force field, this short version of the TM1-containing peptide, PEPT1-29, was observed to exit the membrane, preventing dimerization events. Inspired by a recent study of Martini 3 dimerization of TM-helices [START_REF] Amita R Sahoo | Transmembrane dimers of type 1 receptors sample alternate configurations: Md simulations using coarse grain martini 3 versus alphafold2 multimer[END_REF], we replaced the short sequence with the longer PEPT1-41 variant, which contains additional flanking residues at both the N-and C-termini, effectively reducing the number of ejections. The longer TM1 sequence allowed us to generate 25 successful dimerizations (for each protonation state of Lys716) for subsequent analyses. Because the simulations with Martini 2 and Martini 3 produced quite similar results in terms of dimer associations, and because Martini 3 addresses a number of shortcomings of Martini 2, such as exaggerated protein-protein aggregation, we decided to present only the Martini 3-based results in the remaining sections of this article. The corresponding Martini 2 results can be found in supplementary Figure 9.

Protonation of Lys716 interferes with the formation of a stable TM1-TM2 dimer

The simulations were started with a distance between the two peptides of about 5.2 nm. Most of the time, spontaneous and irreversible dimerization was observed (Figure 1 e1-e2 and f1-f2), but sometimes the peptide containing TM1 was ejected from the membrane (Figure 1 f3-f4-f5). All ejections occurred on the intermembrane side, which corresponds to the lower leaflet in our simulations. When dimerization occurred, we observed that the distance between the two peptides stabilized at 1.5 nm. This event generally occurred within the first 2 µs in each replicate (supplementary Figure 6A). In the Martini 2 simulations, dimerization was irreversible in all runs (supplementary Figure 6A). In contrast, in a considerable number of Martini 3 simulations, TM1 ejection occurred, which affected the dimerization 8/31

process. We therefore had to run up to 63 simulations with Martini 3 to obtain the targeted set of 25 trajectories that ended with a stable dimer (supplementary Figure 6A and6B). For charged Lys716, we observed almost 5 times more ejections (38 ejections versus 25 dimerizations) than for the neutral version (8 ejections versus 25 dimerizations) (see supplementary Figure 6). The lower propensity for dimerization observed with Martini 3 when Lys716 is charged underscores the greater compatibility of a neutral form for self-association within the hydrophobic environment of the membrane.

Protonation of Lys716 acts as a switch for remodeling the TM1-TM2 interface

Next, we analyzed the TM1-TM2 contacts that resulted from the CG simulations after dimerization. To identify the most important dimer conformations from the two sets of 25 trajectories, we clustered the dimers using the GROMOS method [START_REF] Daura | Peptide folding: When simulation meets experiment[END_REF]. This clustering was based on the pairwise RMSD matrix of the backbone beads of each pair of dimerized TM domain conformations. In the Martini 3 simulations with a neutral Lys716, clustering resulted in 5 clusters representing 96% of the conformations, with the first 3 clusters corresponding to 91.4% of the structures. However, when Lys716 was charged, we found that the first cluster was smaller, while the others were larger. Also, the total number of clusters was larger with charged Lys716. This observation is an indication of a greater conformational diversity. The same trend was observed in the Martini 2 simulations (Table 3, Supplementary Figure 9). The interactions between the two TM helices can be observed conveniently by plotting the position of the center of mass of TM1 relative to TM2. All structures used for the clustering analysis were fitted to the backbone beads of TM2, and the positions of the center of mass of TM1 were then plotted in Figure 2A,B. This plot shows the position of TM1 around TM2, which is fixed in the center. Overall, this analysis reveals that the charge state of Lys716 controls the TM1-TM2 association. For charged Lys716, TM1 contacts TM2 on the left and top (Figure 2B), whereas for neutral Lys716, TM1 contacts TM2 on the right and top (Figure 2A). TM-TM contacts are therefore completely different depending on the charge state of Lys716.

The TM domain energy landscape quantifies Lys716-mediated dimer (de)stabilization

A principal component analysis (PCA) was performed on the dimerized structures from the 4 sets of simulations. Figure 2C shows the projection of each conformation onto the first two principal components, in terms of free energy, for neutral Lys716 with Martini 3. For charged Lys716 (as well as for Martini 2 simulations), the PCA plots are shown in supplementary Figure 9. For neutral Lys716, there is one main free energy well (Figure 2C), whereas there are several shallow wells when it is charged (supplementary Figure 9). This observation is consistent with the clustering results (Table 3). The charged residue promotes the exploration of a greater variety of conformations. We compared these results with dimers previously obtained in Ref. [START_REF] Vecchis | A membrane-inserted structural model of the yeast mitofusin fzo1[END_REF] predicted with PREDDIMER [START_REF] Polyansky | Preddimer: a web server for prediction of transmemrane helical dimers[END_REF]. It can be observed that none of the PREDDIMER predictions match the first 5 clusters or the deepest wells (Figure 2A,2B,2C).

For the neutral Lys716 system, the deepest free energy well (Figure 2C) coincides exactly with the center of the first (largest) cluster (Figure 2A). Therefore, we decided to choose this latter structure as our best CG model. We extracted it for further analyses and conformational sampling at the allatom level (see below). This best model shows a compact structure characterized by a crossing angle 9/31 2C) and it is a left-handed antiparallel dimer. The GX 3 G motif contained in TM2 is involved in the interaction with TM1 (Figure 2C), consistent with the literature [START_REF] Brosig | The dimerization motif of the glycophorin a transmembrane segment in membrane: importance of glycine residues[END_REF]. However, this motif is usually observed in right-handed dimers [START_REF] Hubert | Single-spanning transmembrane domains in cell growth and cell-cell interactions: More than meets the eye?[END_REF][START_REF] Walters | Helix packing motifs in membrane proteins[END_REF][START_REF] Zhang | The membrane-and soluble-protein helix-helix interactome: similar geometry via different interactions[END_REF], which is not the preferred arrangement here.

To have a better idea of the contacts, Figure 2D shows a contact-map of this best CG structure, and Figures 2E and2F show a back-mapping of it to an all-atom representation. The part of each helix involved in the interaction ranges from T706 to N720 for TM1 and S738 to L753 for TM2. The N-terminal part of TM1 interacts mainly with the C-terminal part of TM2 (Figure 2B).

Charged Lys716 destabilizes the membrane

The presence of a Lys residue in the middle of a TM helix raises the question of its putative role in the fusion process. Here, we performed CG simulations of the dimerization of TM1 and TM2 (as well as some flanking residues) with charged or neutral Lys716, allowing us to assess the effects of the dimers on the surrounding membrane environment. Thus, we analyzed two parameters associated with the onset of hemifusion, the fusion of the outer layers of each membrane: (i) lipid packing defects, which quantify the hydrophobic membrane surface in contact with the solvent [START_REF] Vamparys | Conical lipids in flat bilayers induce packing defects similar to that induced by positive curvature[END_REF][START_REF] Gautier | Packmem: A versatile tool to compute and visualize interfacial packing defects in lipid bilayers[END_REF] and (ii) lipid tail protrusions, defined as the appearance of a carbon atom (or a coarse-grained bead belonging to a lipid tail) protruding above the level of the phosphate group [START_REF] Larsson | Lipid tail protrusion in simulations predicts fusogenic activity of influenza fusion peptide mutants and conformational models[END_REF][START_REF] Mukarram | Solvent-exposed lipid tail protrusions depend on lipid membrane composition and curvature[END_REF]. When such a protrusion occurs, the hydrophobic lipid tail is assumed to be in contact with the solvent. Thus, both parameters indicate the probability that the hydrophobic tails are exposed to the solvent, which is a prerequisite for hemifusion [START_REF] Pincet | Short-range specific forces are able to induce hemifusion[END_REF][START_REF] Larsson | Lipid tail protrusion in simulations predicts fusogenic activity of influenza fusion peptide mutants and conformational models[END_REF][START_REF] François-Martin | Cooperation of conical and polyunsaturated lipids to regulate initiation and processing of membrane fusion[END_REF].

For lipid packing defects, the π constant provides information on the occurrence and size of such defects, with a higher value indicating more frequent and larger defects. Interestingly, the π constants are higher when the TM domains are present compared to pure membranes (supplementary Figure 7A). The constant for charged Lys716 is slightly higher, which suggests more packing defects, than that for neutral Lys716. The effect is however too modest to draw a definitive conclusion for this parameter. As for the protrusions, the presence of the protein favors their occurrence (Supplementary Figures 7B), which is consistent with the analysis of packing defects. In addition, an increased occurrence of protrusion events is observed when Lys716 is charged, indicating a destabilizing effect of the charge on the membrane. In summary, the route to hemifusion is favored not only by the presence of the TM domain, but the charge state of Lys716 also matters.

Replica exchange MD provides a refined atomistic model and verifies stable TM assembly

To fine-tune a model of the whole TM domain, we performed all-atom simulations (AA) using an enhanced sampling technique, temperature replica-exchange molecular dynamics (T-REMD). As a starting structure, we chose the best dimer of the Martini 3 simulations obtained with neutral Lys716, as described above and in Figure 2. The missing loop was added and sampled as described in the supplementary methods. Finally, the system was back-mapped to an AA representation (see Figure 2E and2F) and used as the initial structure for a T-REMD simulation of 500 ns. In total, 38 replica were simulated ranging from 303.15 K to 400 K. This REMD protocol allowed us to test the robustness of the model by exposing it to high temperatures while maintaining a physical Boltzmann distribution at each temperature.

We first analyzed the behavior of each replica as it propagated through temperature space. Of the 38 replicas, 19 showed persistent interactions between the two TM helices (center of mass / center of mass distance below 2 nm) throughout the simulation, 12 showed TM dissociations, while 7 showed TM dissociations and re-associations (Supplementary Figure 10A). Importantly, the Cα RMSD of TM1 and TM2 showed minimal structural changes for each replica when the two TM regions maintained their interaction during the simulation (Supplementary Figure 10B), indicating stability and robustness of the original structure.

Next, we focused on the ensemble of structures at the lowest temperature (303.15 K) (Figure 3). In total, 79.2 % of the conformations presented stable TM1 / TM2 contacts. In addition, all conformations showed stability of their secondary structures in the two TM helices, even those which displayed 11/31 dissociation (Supplementary Figure 11). Figure 3A shows the position of TM1 around TM2 fixed at the center, similar to Figure 2A for CG simulations, in terms of free energy. The conformations of the bottom temperature quickly left the area of the starting structure (represented by a cross), and populated mainly the right side of TM2, as well as its upper side but in a less pronounced way. The Martini 3 simulations had captured the right and upper parts overall, but not the exact upper right region found here. [START_REF] Vecchis | A membrane-inserted structural model of the yeast mitofusin fzo1[END_REF]. TM1 and TM2 are seen in yellow, the other residues used for the REMD are seen in orange, HRN is in purple, HR1 in blue, HR2 in green and the domain GTPase is in red. The residues that do not belong to a specific domain are in light grey. The second image is the TM portion, zoomed in. The residues before Ser702 and after Lys761 are represented in transparent gray.

In the following, we consider only the conformations with TM1 / TM2 in contact. From this ensemble 12/31 of structures, we extracted the two best models of the entire TM domain from the two wells with the lowest free energy, as shown in Figure 3A-F. The two wells are labeled C and D in panels 3A and 3B, whose structures are shown in panels 3C and 3D. In both models C and D, TM1 touches TM2 in the lower part (the most likely contacts in Figure 3E-F are below the diagonal of the contact map), especially in model C. Both models C and D were then compared with the overall Fzo1 model that we proposed in 2017 [START_REF] Vecchis | A membrane-inserted structural model of the yeast mitofusin fzo1[END_REF] to check which model fits better (Figure 3G-H). As can be seen, model D (Figure 3D/H) fits better with the overall model of Fzo1 because the TM domain can be connected to the rest of the model at both the N-and C-terminal sides. We therefore selected model D as our best all-atom model. However, it is important to keep in mind that the TM contacts have a certain degree of flexibility, so that the alternative conformation C cannot be completely excluded.

We also analysed the promiscuous stabilizing interactions between TM1 and TM2 (Figure 3E-F). While we observed a similar set of residues involved in the interactions in TM1, there is a clear shift in the interactions involving TM2 (Figure 3). Instead of Gly745 and Gly749 of the GX 3 G motif, residues Val741 and Leu748 were the most involved in the interactions. We also found that among the structures located within the energy wells (Figure 3), Thr706, Leu709, Gly713, Lys716 from TM1 and Leu737, Val741, Val744, Leu748 from TM2 were the most promiscuous residues.

In terms of hydrogen bonds, we found some between the side chains of Lys716 and Ser738 with 8% of persistence. However, the latter residue frequently formed hydrogen bonds with the side chain of the TM1 residue Asn720 (12% of persistence) as well as with the residues of the loop from Lys726 to Lys735 (both side chains and backbone with 7% of persistence) and most frequently with the solvent (21% of persistence).

Experiments show that the polar residue 716 is crucial for mitochondrial respiration in yeast

We sought to obtain experimental confirmation of our theoretic observations on residue K716. We reasoned that mutating this Lysine into apolar residues (I or V) should impact mitochondrial fusion efficiency while mutation into polar (S) or charged (R) residues should have no effect. To test this prediction, we took advantage of the established link between mitochondrial fusion efficiency and respiratory growth of Saccharomyces cerevisiae. Depending on the carbon source provided in their media, yeast cells grow through either fermentation or respiration. Dextrose is a fermentative carbon source that inhibits respiration whereas Glycerol is a fully respiratory carbon entry [START_REF] Rolland | Glucose-sensing and-signalling mechanisms in yeast[END_REF][START_REF] Gancedo | Yeast carbon catabolite repression[END_REF]. Consequently, since mitochondrial fusion is essential for respiration, its inhibition abolishes yeast growth in Glycerol-containing media [START_REF] Hermann | Mitochondrial fusion in yeast requires the transmembrane gtpase fzo1p[END_REF]. WT, fzo1∆ and cells expressing distinct versions of Fzo1 mutated in K716 as the sole source of Fzo1 (Fig. 4a) were thus subjected to serial dilutions growth assays on Dextrose or Glycerol-containing media at 23, 30 or 37°C (Fig. 4b). All mutants were expressed at levels comparable to WT Fzo1 (Fig. 4A). However, respiratory growth of FZO1 K716I and FZO1 K716V cells was abolished at all temperatures, similar to fzo1∆ cells. In contrast, respiratory growth of FZO1 K716S and FZO1 K716R cells was not affected, similar to WT cells. These results indicate that mutation of Lys716 to an apolar residue (K716I or K716V) impairs respiration likely because of an inhibition of mitochondrial fusion. Moreover, mutation of Lys716 to an Arg (K716R), the main alternative in Mammalian mitofusins, or Mutation of Lys716 to a polar residue (K716S) has no effect on respiration. In this context, the positive charge of Arg or Lys does not seem to be mandatory, but the residue polarity of Lys, Arg or Ser at position 716 would be essential for mitochondrial fusion.

Discussion

Our extensive simulations provide new insights into the structure and dynamics of the transmembrane domain of Fzo1. Comparison with previous predictions reveals important differences, while independent validation suggests that the refined model accurately reproduces the folding of the TM domain. In this discussion, we first examine the robustness of the predicted TM1-TM2 interface in different simulation methods. Comparison of coarse-grained and all-atom approaches shows the robustness of the main TM 13/31 An independent prediction of the Fzo1-Ugo1 complex from AlphaFold2 is used to externally validate our new TM domain fold. We examine the remarkable agreement between these completely different approaches. We then discuss the functional implications of our refined structure, including how specific residues may be actively involved in the membrane fusion process. Our model sheds light on how the properties of the TM domain may facilitate the hemifusion and full fusion stages. Finally, we compare the results with other studies of dimerization, which give a consistent picture of such processes when examined with coarse-grained MD simulations.

The TM1-TM2 interface is robust across force fields and levels of representation

To investigate the possible association of TM1 and TM2, we generated twenty-five 10 µs coarse-grained (CG) MD trajectories showing dimerization. Although the TM1-TM2 biological construct is connected by a loop in the intermembrane space, we intentionally cut the TM domain into two peptides so that 14/31 dimerization was freely driven by TM1-TM2 contacts within the membrane. This assembly approach can be considered equivalent to docking, but confined within an explicit membrane environment. Because we used the CG Martini force fields [START_REF] Siewert J Marrink | The martini force field: coarse grained model for biomolecular simulations[END_REF][START_REF] Monticelli | The martini coarse-grained force field: extension to proteins[END_REF][START_REF] Paulo | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] with smoother energy surfaces than a fully atomistic model, we were able to benefit from accelerated kinetics and thus efficiently sample the different modes of interaction of the two TM helices. As a result, we obtained robust statistics on TM1-TM2 contacts.

protocol was first tested with the Martini 2 force field (supplementary Figure 9). During the progress of this work the Martini 3 force field was released, which encouraged us to use it as well. Interestingly, the Martini 2 simulations yielded a very similar dimerization pattern compared to Martini 3, albeit with a stronger association, which is a known artefact of Martini 2 [START_REF] Paulo | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF]. Conformational clustering showed the same trends for both force fields, i.e., the charged Lys716 system has a larger number of clusters, and the first cluster of the simulations with neutral Lys716 accounts for more than 50% of the sampled structures. Experimentally, some 2 H experiments were performed on a TM helix containing a Lys in its center [START_REF] Gleason | Buried lysine, but not arginine, titrates and alters transmembrane helix tilt[END_REF]. A charged Lys was leading to a multistate behavior in terms of peptide orientation in the membrane (tilt and azimuthal rotation), whereas a neutral Lys gave a single orientation. In our CG simulations we have two TM helices dimerized, but we find the same kind of behavior since we get a single free energy well for neutral Lys716, and multiple shallow ones for the charged version.

When the model was subjected to all-atom simulations (AA) using an enhanced sampling method (T-REMD), we observed a slight shift in the most frequent interactions between the two TM helices. Lamprakis et al. showed that Martini 3 does not always favor the experimentally solved interface of interactions between TM domains [START_REF] Lamprakis | Evaluating the efficiency of the martini force field to study protein dimerization in aqueous and membrane environments[END_REF] and recommended that a refinement procedure be used. Nevertheless, the interfaces of the interactions involving TM1 are overall the same in the Martini 3 and CHARMM36 refinements. The change to a detailed AA representation affected mainly the interactions involving TM2, while in the overall model, the interface of the interactions and the crossing angle were only slightly affected.

How can we interpret TM1 expulsion from the membrane with Martini 3?

Next to experiments of biophysics or cellular biology, CG simulations with Martini 2 or 3 have been widely used to predict TM-TM helix dimerization (or higher order oligomerization), for example in recent works [START_REF] Ostuni | Cx3cl1 homo-oligomerization drives cell-to-cell adherence[END_REF][START_REF] Jackson | The guidance and adhesion protein flrt2 dimerizes in cis via dual small-x3-small transmembrane motifs[END_REF][START_REF] Alavizargar | Lipidmediated association of the slg1 transmembrane domains in yeast plasma membranes[END_REF][START_REF] Valério | Parainfluenza fusion peptide promotes membrane fusion by assembling into oligomeric porelike structures[END_REF]. Both force field versions are able to predict TM-TM interfaces, but Martini 2 has a tendency to overaggregate TM segments. Considerable efforts have thus been put in the development of Martini 3 to solve this issue [START_REF] Paulo | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF]. This new version was recently tested on many known TM-TM dimers with success [START_REF] Amita R Sahoo | Transmembrane dimers of type 1 receptors sample alternate configurations: Md simulations using coarse grain martini 3 versus alphafold2 multimer[END_REF], but the authors highlighted the ability of Martini 3 to sample alternative conformations as well as the importance of flanking residues (around TM helices). Subsequent studies found the possible ejection of TM helices with Martini 3 [START_REF] Valério | Parainfluenza fusion peptide promotes membrane fusion by assembling into oligomeric porelike structures[END_REF][START_REF] Spinti | Room for improvement in the initial martini 3 parameterization of peptide interactions[END_REF], ending adsorbed on the membrane in a horizontal orientation. To avoid this phenomenon, some solutions were proposed such as adding position restraints or rescaling protein-water van der Waals interactions [START_REF] Claveras Cabezudo | Scaling protein-water interactions in the martini 3 coarse-grained force field to simulate transmembrane helix dimers in different lipid environments[END_REF]. In summary, Martini 3 is still in a testing phase by the scientific community. We are learning progressively how it behaves and how to interpret its outcomes by testing it on several different systems. In our case, we tested the rescaling procedure [START_REF] Claveras Cabezudo | Scaling protein-water interactions in the martini 3 coarse-grained force field to simulate transmembrane helix dimers in different lipid environments[END_REF], but still observed TM1 ejections (data not shown).

In this context, how can we interpret the ejections of TM1 observed in our simulations with Martini 3? On one hand, this is unexpected since we are dealing with a TM helix. 2 H NMR experiments have confirmed that a single TM helix containing a Lys in its center assumes a transmembrane topology regardless of the charge state [START_REF] Gleason | Buried lysine, but not arginine, titrates and alters transmembrane helix tilt[END_REF]. On the other hand, we have to consider that there are 3 consecutive polar residues STK in the middle of TM1. Moreover, Lys 716 may be charged depending on the conditions. In the context of the full protein, which bears the HR1 cytoplasmic domain on its N-terminal side, these ejections are probably not realistic. But due to the simplified description of coarse-grained vs all-atom representation (fixed secondary structure, etc.), due to the absence of the cytoplasmic domains, this is how Martini 3 is warning us that something is going on with TM1: it is not an ideal helix for a hydrophobic environment like the center of a membrane. Yet, we demonstrated that TM1 destabilizes the membrane. If Lys716, were to be charged, this would be even more pronounced. This is interesting, because in the context of membrane fusion, a TM helix able to destabilize the membrane is clearly 15/31 an asset, such as in SNARE proteins [START_REF] Jelger Risselada | Caught in the act: Visualization of snare-mediated fusion events in molecular detail[END_REF][START_REF] Han | Synaptobrevin transmembrane domain determines the structure and dynamics of the snare motif and the linker region[END_REF]. We discuss this aspect further 2 sections below. Last, as stated by Sahoo et al., we observed that the flanking residues matter, since the longer PEPT1-41 underwent less ejections than the shorter PEPT1-29 [START_REF] Amita R Sahoo | Transmembrane dimers of type 1 receptors sample alternate configurations: Md simulations using coarse grain martini 3 versus alphafold2 multimer[END_REF]. We also observed that charged Lys716 was ending in far more ejections than neutral Lys716. The ratio of ejections may thus be seen as a proxy towards the likeliness to the membrane.

Extensive multiscale simulations yield a new structural model of the Fzo1 TM domain improving over previous predictions

The TM domain model refined with our multiscale protocol shows significant deviations from the 2017 prediction [START_REF] Vecchis | A membrane-inserted structural model of the yeast mitofusin fzo1[END_REF]. The soluble part of the previous model could be validated by experimental mutation studies, while the TM1-TM2 contacts of the transmembrane part were based on an ab initio prediction by the PREDDIMER method [START_REF] Polyansky | Preddimer: a web server for prediction of transmemrane helical dimers[END_REF] that could not be experimentally assessed at that time. Of the three PREDDIMER models created, the highest scoring conformation was selected, a right-handed antiparallel dimer characterized by a crossing angle of 119.7 • . In contrast, the model created here with Martini 3 is a left-handed antiparallel dimer with a crossing angle of -137.4 • , and the model refined with REMD has an angle of -161.4 • . In addition, the Martini 3 model shows involvement of the GX 3 G motif, while the refined model no longer shows interaction of these residues with TM1. This motif could thereby be free for other contacts such as homodimerization, or with another partner in the membrane. In addition, TM1 shows no matching interactions with the residues observed with PREDDIMER. T715 is most closely associated with the GX 3 G motif in the PREDDIMER model, whereas in the new model it is not T715 but L709, Y710, and G713, among others.

An independent Fzo1-Ugo1 complex prediction matches the novel TM domain fold

The field of structure prediction has been revolutionized recently by the success of AlphaFold2 (AF2) [START_REF] Jumper | Highly accurate protein structure prediction with alphafold[END_REF][START_REF] Evans | Protein complex prediction with alphafold-multimer[END_REF]. This observation prompted us to use AF2 to independently predict the structure of the TM domain of Fzo1.

Using AF2, we first predicted the TM domain of Fzo1 alone, then a whole Fzo1 monomer and a homodimer of Fzo1. All these predictions led to the two TM helices unstructured, with a low pLDDT (Predicted Local Distance Difference Test), which is the confidence value per residue (Supplementary Figure 13). This observations seems to echo the difficulty to predict membrane-inserted protein segments with AF2.

In yeast, Fzo1 has an important biological partner in the outer membrane named Ugo1 [START_REF] Sesaki | UGO1 Encodes an Outer Membrane Protein Required for Mitochondrial Fusion[END_REF] which is involved in mitochondrial fusion [100, 101]. Interestingly, when two Fzo1 monomers and two Ugo1 monomers were subjected to AF2, we observed that the two TM helices were correctly predicted to interact with each other (Figure 5). Strikingly, compared to our best model from the REMD simulations (Figure 3D), the two structures are very close: the RMSD (on backbone Cα atoms of TM1 and TM2 only) is as low as 0.27 nm, the promiscuous residues are very similar, and the crossing angles are close (-159 • for AF2, -163 • for our model) (Figure 5). Moreover, the AF2 model shows that the TM domain of each Fzo1 monomer interacts with each other around the TM helices (mainly TM2-TM2 contacts). This result is interesting because it suggests that Fzo1-Fzo1 dimerization may also involve contacts between the two TM domains. As the GX3G motif does not actually interact with TM1, these TM2 residues could be involved with the other monomer, even if we find closer interactions with other residues.

In summary, the overall agreement between our strategy and an artificial-intelligence-based method reinforces the pertinence of our physics-based model.

Could the Fzo1 TM domain play an active role in membrane fusion?

The fusion process starts with the two membranes getting close to one another. The approach is followed by the formation of a stalk intermediate resulting from the mixing of the outer bilayer leaflets. Then 16/31 Recent experimental studies have shown the importance of lipid conformations in the initiation of the outer leaflets mixing. The ability of lipids to splay, i.e. expose a part of their hydrophobic tail towards the solvent also known as protrusion, was shown to trigger membrane fusion [103]. It was also shown that lipid packing defects, small hydrophobic areas exposed to the solvent, were qualitatively correlated to the nucleation rate of fusion [START_REF] François-Martin | Cooperation of conical and polyunsaturated lipids to regulate initiation and processing of membrane fusion[END_REF]. Protrusions and packing defects are logically directly correlated on lipid composition and membrane curvature [104,105,[START_REF] Gautier | Packmem: A versatile tool to compute and visualize interfacial packing defects in lipid bilayers[END_REF][START_REF] Mukarram | Solvent-exposed lipid tail protrusions depend on lipid membrane composition and curvature[END_REF]. These two parameters are thus critical in the initiation of membrane fusion. membrane fusion mediated by proteins, it is now well established that their TM domains play a pivotal role [106]. In the well-studied case of SNARE proteins, the unique presence of a single TM helix in each leaflet catalyses outer leaflet mixing [START_REF] Jelger Risselada | Caught in the act: Visualization of snare-mediated fusion events in molecular detail[END_REF]. The TM helix sequence matters since some mutants induce less protrusions than the wild type [START_REF] Han | Synaptobrevin transmembrane domain determines the structure and dynamics of the snare motif and the linker region[END_REF]. In addition, the lipids are more perturbed near the TM helix than far away. Another well known example is the influenza fusion peptide [START_REF] Larsson | Lipid tail protrusion in simulations predicts fusogenic activity of influenza fusion peptide mutants and conformational models[END_REF]. Again, the TM sequence is decisive for the generation of protrusions around it (leading to membrane fusion), while some mutants are less efficient. A recent experimental study on model TM helices has also shown that the amount of lipid splay (i.e. protrusion) was directly correlated with membrane fusion [103]. In this example, membrane fusion also depended on the amino-acid sequence, a poly LV16 being more efficient than a poly L16 TM peptide. Similar to these examples, our results suggest a destabilization of the membrane promoted by the presence of the Fzo1 TM domain.

We discuss now the possible role of Lys716 in membrane fusion. First, this lysine is well conserved in fungi (Figure 13). In mammals, an arginine is rather found at this position. In organisms that have a full TM domain (not a paddle which only partially binds the membrane like in bacterial BDLP), it is interesting to note that evolution have conserved a basic residue in the middle of TM1, which raises questions about its role. Our experiments showed that mitochondrial respiration is disabled when Lys716 is mutated to an apolar residue (Leu or Val), but mutation to an Arg or Ser has no effect on respiration. Amino-acid polarity is thus required, but having a basic one is not mandatory. These data are in line with the membrane destabilization hypothesis, since a polar residue will anyway pertub the membrane (recall there are two other polar residues next to Lys716). However, these findings do not preclude the possibility of Lys716 positive charge to play a role. The experiments of Gleason et al. evaluated the pKa of a Lys in the middle of a TM helix to be around 6.2 at 50 • C [START_REF] Gleason | Buried lysine, but not arginine, titrates and alters transmembrane helix tilt[END_REF]. Correcting for temperature, the pKa was then estimated to be ∼ 6.5 at 37 • C and ∼ 6.8 at 25 • C. Given an acidic pH in yeast cytoplasm[107], these values would suggest that Lys716 could be, at least partially, protonated. The role of positively charged residues is not new in the context of fusion. In the case of the SNARE fusion machinery, Lindau et al. examined the transmembrane domains in detail and hypothesized that well-placed charges could significantly assist membrane deformation and destructuring, thereby initiating the fusion process. In particular, they investigated the role of charged residues at the end of a Syb2 TM domain construct when pulled into the membrane [START_REF] Lindau | Coarse-Grain Simulations Reveal Movement of the Synaptobrevin C-Terminus in Response to Piconewton Forces[END_REF]. It has further been suggested that charged motifs play an important role for fusion to occur, such as shown by the bilayer destabilization by a conserved membrane-embedded motif at the juxtamembrane region of the vesicle-anchored v-SNARE comprising several basic residues [108]. In TM1 of Fzo1, the Lys is in the middle of the TM helix, not flanking it, but it does not exclude that a similar mechanism may be operating. Overall, these observations echo the idea of a protonation-state induced switch of Lys716 to assist membrane fusion. However, the lack of effect when the Lys is mutated to a polar, but not titratable residue, is arguing against a direct effect on mitochondria fusion. Therefore, the ability of Lys716 to be protonated could be associated with the fusion of other membranes to which Fzo1 is associated. Besides mitochondrial fusion, we cannot exclude that this titration could be associated to another function.

Conclusion

We have constructed an improved model of the TM domain of mitofusins using both coarse-grained and all-atom molecular dynamics. This model has been further confirmed using the deep-learning tool AlphaFold2. This model has revealed the role of Lys716 that is found in the middle of the membrane. The importance of this residue is confirmed by its evolutionary conservation and the effect of its mutations on mitochondrial fusion. Further studies will therefore be needed to explore the exact role of this residue in mitochondrial fusion but could also target other functions associated with mitofusins.

Figure 1 .

 1 Figure 1. Protocol. (A) Sequence definition of the peptides PEPT1-29, PEPT1-41 and PEPT2. The TM domains are represented in yellow and the flanking residues are in orange. This color code is maintained for the snapshots in B,C,E,F. (B) Images of the all atom structures of the peptides PEPT1-41 and PEPT2. (C) Snapshot of the coarse-grained system. In blue are the phosphate groups, in silver are the lipids POPC and POPE. (D) A detailed scheme of the protocol. (E) Center of mass distances throughout the simulation, with snapshots of the beginning and end point. (F) Center of mass distances throughout the simulation, for a replica with TM1 expelled from the membrane. Snapshots throughout the simulations.

  fzo1∆::LEU2 pRS416-FZO1 fzo1∆ (MCY620) MATα ura3-1 trp1-1 leu2-3,112 his3-11,15 can1-100 This study fzo1∆::LEU2 pRS314 FZO1 (MCY1261) MATα ura3-1 trp1-1 leu2-3,112 his3-11,15 can1-100 This study fzo1∆::LEU2 pRS314-FZO1 FZO1-K716R (MCY1402) MATα ura3-1 trp1-1 leu2-3,112 his3-11,15 can1-100 This study fzo1∆::LEU2 pRS314-FZO1-K716R FZO1-K716S (MCY1425) MATα ura3-1 trp1-1 leu2-3,112 his3-11,15 can1-100 This study fzo1∆::LEU2 pRS314-FZO1-K716S FZO1-K716V (MCY2258) MATα ura3-1 trp1-1 leu2-3,112 his3-11,15 can1-100 This study fzo1∆::LEU2 pRS314-FZO1-K716V FZO1-K716I (MCY2259) MATα ura3-1 trp1-1 leu2-3,112 his3-11,15 can1-100 This study fzo1∆::LEU2 pRS314-FZO1-K716I

Figure 2 .

 2 Figure 2. TM1-TM2 contacts determined using Martini 3 and best model extraction. (A,B) Positions of the center of mass of TM1 around the center of mass of TM2 (midpoint) for neutral (A) and charged (B) Lys716. The coloring of TM1 positions is based on the cluster to which the structure belongs to. Shown on the right side of each plot are the central structures of the first three clusters, colored similarly. (C) Free energy projection on the first two principal components of a PCA analysis (neutral Lys716 system). This PCA was performed on a concatenated trajectory containing all CG simulations (Martini 2 and 3, neutral and charged Lys716). The crosses 1, 2 and 3 correspond to the 3 dimers obtained with PREDDIMER in ref. [30]. Shown on the right are two snapshots of the best model (labeled "a" and corresponding to the center of cluster 1 in panel A and to the deepest free energy well), together with its crossing angle. The glycines of the GX 3 G motif are shown in van der Waals representation. Lys716 is shown in stick representation. (D) Contact map of the best model. (E,F) Best model back-mapped to an all-atom representation with the loop between TM1 and TM2 added. The protein is shown in a cartoon representation. TM1 and TM2 are shown in yellow and the flanking residues (with the loop) are shown in orange. Lys716 is shown in stick representation, and the GX 3 G motifs are shown in van der Waals representation. The phosphorus atoms of the lipids POPC and POPE are shown as blue spheres, and the water molecules are shown as cyan sticks. The loop was added according to the procedure described in the supplementary material.

Figure 3 .

 3 Figure 3. TM domain prediction based on REMD simulations. (A) Positions of the center of mass of TM1 around the center of mass of TM2 (center point) colored as function of the free energy. The orange cross is the best CG model obtained with Martini 3 (in Figure 2) used as a starting point of the REMD. (B) Free energy projection on the first two principal components of a PCA analysis. This PCA was performed on a concatenated trajectory containing Martini 3 (neutral K716) and all-atom (bottom temperature) conformations. On the plot, only the all-atom data are shown. The crosses 1, 2 and 3 correspond to the 3 dimers obtained with PREDDIMER in ref. [30]. The structures selected for both (A, B) have a TM1-TM2 COM-COM (center-of-mass) distance of less than 2 nm. (C,D) Center structure of the first and third cluster respectively, for a cutoff of 0.2 nm. These clusters roughly correspond to the 2 main free energy wells (bluest parts) seen in (A,B) and labeled C and D. The TM domains and the flanking residues are in cartoon representation. TM1 and TM2 are in yellow, and the flanking residues (with the loop) are in orange. The glycines of the GX 3 G motif are shown in van der Waals representation. The Lys716 is shown in stick representation. (E,F) Contact map of the models (C,D) respectively. (G,H) Addition of respectively (C) and (D) to the model of Fzo1 built by De Vecchis et al. in 2017[START_REF] Vecchis | A membrane-inserted structural model of the yeast mitofusin fzo1[END_REF]. TM1 and TM2 are seen in yellow, the other residues used for the REMD are seen in orange, HRN is in purple, HR1 in blue, HR2 in green and the domain GTPase is in red. The residues that do not belong to a specific domain are in light grey. The second image is the TM portion, zoomed in. The residues before Ser702 and after Lys761 are represented in transparent gray.

Figure 4 .

 4 Figure 4. The polarity of K716 is essential for Fzo1 function. A) Anti-Fzo1 and anti-Pgk1 immunoblots of total protein extracts prepared from, fzo1∆ shuffle strains transformed with an empty plasmid or plasmids expressing FZO1 WT or mutated in K716. B) Dextrose and glycerol growth spot assay with strains used in A

Figure 5 .

 5 Figure 5. AlphaFold2 prediction of Fzo1 in interaction with Ugo1 (A) Scheme of the protocol. (B,C,D,E) The four proteins are in cartoon representation. In dark blue and light blue are the models of Ugo1. TM1 and TM2 are seen in yellow, the other residues used for the REMD are seen in orange, HR1 in blue and HR2 in green. The residues that do not belong to a specific domain are in light grey. (B) The overall of AlphaFold. (C) Zoom on the TM domains of the Fzo1-small dimers. (D) One Fzo1-small monomer. (E) Zoom on TM1 and TM2 of Fzo1-small monomer. (F) Contact map of TM1 and TM2.

Figure 6 .

 6 Figure 6. Distance between the center of mass of the two TM domains, as a function of time. (A) Results for Martini 2 and Martini 3 systems. (B) Results for Martini 3 simulations, where TM1 exiting the membrane was observed.

Figure 7 .

 7 Figure 7. Caption

Figure 8 .

 8 Figure 8. Position of TM1 relative to TM2 colored as a function of the free energy. (A) Lys716 charged. (B) Lys716 neutral.

Figure 9 .

 9 Figure 9. PCA colored as a function of the free energy.

Figure 10 .

 10 Figure 10. Supplementary results of the REMD. (A) Distance between the center of mass of TM1 and TM2 as a function of time for for all replicas. (B) RMSD of TM1 and TM2 (considering C α only) with respect to the starting structure as a function of time.

Figure 11 .

 11 Figure 11. Supplementary results of the REMD. Secondary structure of the ensemble of structures at 303.15 K.

Figure 12 .

 12 Figure 12. PCA results of the REMD. (A) Martini 3 Lys716 neutral structures colored as function of free energy. (B) REMD structures colored as function of free energy. (C) REMD structures colored as a function of the cluster number to which they belong to.

Figure 13 .

 13 Figure 13. Fzo1 model of Alphafold. The model is shown in cartoon representation. (A) Colored as a function of plddt. (B) The domains are colored, and the residues outside of the domains are in grey. (C) Model built in 2017 [30], for comparison purposes.

Table 1 .

 1 Plasmids used in this study.

	Name (Collection number)	Description	Reference
	pRS314 (MC219)	CEN, TRP1, Amp	[79]
	pRS314-FZO1 (MC250)	CEN, FZO1 promoter-FZO1, TRP1, Amp	[80]
	pRS314-FZO1-K716R (MC430) CEN, FZO1 promoter-FZO1-K716R, TRP1, Amp	This study
	pRS314-FZO1-K716S (MC437) CEN, FZO1 promoter-FZO1-K716S, TRP1, Amp	This study
	pRS314-FZO1-K716V (MC587) CEN, FZO1 promoter-FZO1-K716V, TRP1, Amp This study
	pRS314-FZO1-K716I (MC589)	CEN, FZO1 promoter-FZO1-K716I, TRP1, Amp	This study

Table 2 .

 2 Saccharomyces cerevisiae strains used in this study.

	Name	Genotype	Reference
	FZO1 (MCY571)	MATα ura3-1 trp1-1 leu2-3,112 his3-11,15 can1-100	[30]

Table 3 .

 3 Cluster populations from the CG Martini simulations. Shown is the fraction of the 5 first clusters, the sum of these 5 fractions and the total number of clusters.

	FF ver-sion	charged state of Lys716	cluster 1 (%)	cluster 2 (%)	cluster 3 (%)	cluster 4 (%)	cluster 5 (%)	total (%)	total number of clusters
	Martini 3 neutral	55.7	21.6	14.1	3.4	1.9	96.7	16
	Martini 3 charged	34.1	26.7	19.3	6.7	2.3	89.1	25
	Martini 2 neutral	60.0	13.9	10.0	7.4	3.0	94.3	13
	Martini 2 charged	31.9	25.7	20.6	10.1	4.3	92.6	17

Acknowledgments

This work was supported by Agence Nationale de la Recherche ANR-19-CE11-0018, MITOFUSION and by the "Initiative d'Excellence" program from the French State (Grant "DYNAMO", ANR-11-LABX-0011). This work was granted access to the HPC resources of IDRIS under the allocation 2022-A0130701714 made by GENCI to MB. Computational work was performed using HPC resources from LBT-HPC thanks to support from Geoffrey Letessier. Drs Cathy Etchebest, Nicolas Floquet and Miguel Machuqueiro are acknowledged for stimulating discussions. Lucas Guémené is thanked for his help with experiments and Hanene Samia Ameur for testing Martini3 rescaling schemes.

Data availability

R-markdown notebook and Zenodo archive.

7/31

Author Contributions

AT and PFJF contributed equally to this work. RV, AT and PFJF designed the research. RV carried out all simulations. RV, MB, AT and PFJF analyzed the simulation data. LC carried out the experiments. LC and MMC analysed the experiments. RV, MB, AT and PFJF wrote the article.

Supplementary Material

Coarse-grained loop conformational sampling

System Building

The center of the first cluster obtained through Martini 3 simulations with a neutral Lys716 was selected for the following simulations. The lipids and solvent of the corresponding frame were also extracted. This dimer conformation was used as a starting point for the addition of the missing residues between PEPT1-41 and PEPT2 (containing respectively TM1 and TM2).

In order to add the missing residues between TM1 and TM2, a structure of the loop was built using Basic Builder 1.0.2 on the Mobyle server (https://mobyle.rpbs.univ-paris-diderot.fr) [START_REF] Néron | Mobyle: a new full web bioinformatics framework[END_REF]. This structure was then converted to coarse-grained (CG) representation using the martinize2.py workflow module of the Martini 3 force field. As a first step, the residues at the end of the loop were used to align the structures of the dimer with the structure of the loop. Coordinates of residues Lys727 to Lys736 of the dimers were then removed and replaced by the coordinates of the structure available. Residues Arg695 to Thr701 included were also removed.

Simulation parameters

The GROMACS 2018.5 program was used to perform all MD simulations [START_REF] James | Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF]. This system was then subjected to two energy minimization phase of 5000 and 10000 steps. Position restraints of 100000 kJ.mol -1 .nm -2 were used for the TM1 and TM2 residues. Subsequently, all simulations were performed within the NPT ensemble. The velocity-rescaling thermostat [START_REF] Bussia | Canonical sampling through velocity rescaling[END_REF] at 303.5K and the Berendsen barostat at 1bar [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] was applied for 500 ps of equilibration, with a time step of 0.01 ps. The position restraints for the TM1 and TM2 residues were lowered at 20000 kJ.mol -1 .nm -2 . A production run of 5 $µ$s followed at 303.5K using the velocity-rescaling thermostat [START_REF] Bussia | Canonical sampling through velocity rescaling[END_REF] (lipid, water and proteins coupled separately) and at 1 bar using the Parrinello-Rahman barostat [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] (compressibility of 3.0 x 10\ -4 bar). Pressure coupling was applied semi-isotropically. The position restraints for the same residues were fixed at 1000 kJ.mol -1 .nm -2 . A time step of 0.01 ps was used with the leapfrog integrator. Lennard-Jones interactions were cutoff at 1.1 nm. Bond lengths were constrained using the LINCS algorithm [START_REF] Hess | LINCS: A linear constraint solver for molecular simulations[END_REF]. The reaction-field method [START_REF] Ilario | A generalized reaction field method for molecular dynamics simulations[END_REF] was used for evaluating electrostatic interactions, with a Coulomb distance cutoff of 1.1 nm, a relative dielectric constant of 15. The neighbor list was updated every 20 steps.

At the end of this stage, we thus have a coarse-grained system with TM1 and TM2 connected by the loop and embedded in the bilayer surounded by solvent and ions.

Clustering

The loop was submitted to a conformation based clustering. The first 200 ns of the production were ignored, resulting in a total of 9601 frames. The frames underwent a clustering with the GROMOS method [START_REF] Daura | Peptide folding: When simulation meets experiment[END_REF] using a cutoff of 0.2 nm (in order to obtain around 5 clusters representing 80-90 % of all conformations). The center conformation of the first cluster was then selected for the next step.

Supplementary Results

26/31