Structural analysis of M1AP variants associated with severely impaired spermatogenesis causing male infertility

Evolutionary conservations

It is known that structurally and functionally important amino acids are mostly conserved during evolution (Ashkenazy et al., 2016). To reveal the importance of variant regions, we aligned the UniProt (The UniProt Consortium, 2019) protein sequences of M1AP orthologs from human (Homo sapiens, UniProt ID: Q8TC57), chimpanzee (Pan troglodytes, UniProt ID: H2R3U8), cat (Felis catus, UniProt ID: A0A337RXR8), wild boar (Sus scrofa, UniProt ID: A0A5G2QVF2), bovine (Bos taurus, UniProt ID: E1BF42) and house mouse (Mus musculus, UniProt ID: Q9Z0E1) by using T-Coffee with its default parameters (Notredame, Higgins & Heringa, 2000; Di Tommaso et al., 2011). The alignment was visualized via ESPript (v3.0) (Robert & Gouet, 2014).

Molecular modeling and validation

Sequence of the canonical isoform of M1AP protein (UniProt ID: Q8TC57-1) was used to build models. GalaxyWeb (Ko et al., 2012), I-TASSER (Roy, Kucukural & Zhang, 2010), Phyre2 (Kelley et al., 2015), PRIMO (Hatherley et al., 2016), RaptorX (Peng & Xu, 2011) and Robetta (Kim, Chivian & Baker, 2004) were used to model different structures with their default parameters. Regular comparative algorithm (we called Robetta model) and domain prediction algorithm (we called Robetta-domain model), which are part of Robetta webserver (https://robetta.bakerlab.org), were applied separately (Kim, Chivian & Baker, 2004). Robetta algorithms also applied an ab initio prediction for part of the structure and compared their reliability with its template-based predictions.

A list of the templates automatically selected by the modeling algorithms, along with their experimental methodologies, resolution values, and Ramachandran outliers are available in Table S1. Since all structure predictors do not list the identity, similarity, and coverage details, we aligned the templates with M1AP via NCBI’s Protein BLAST (https://blast.ncbi.nlm.nih.gov/) by applying default parameters except the removed E-value threshold (Table S2).

Quality of the structures were assessed by Verify3D (Eisenberg, Lüthy & Bowie, 1997), ERRAT (Colovos & Yeates, 1993), Prove (Pontius, Richelle & Wodak, 1996), PROCHECK (Laskowski et al., 1993) and WHATCHECK (Hooft et al., 1996) algorithms on SAVES v6.0 webserver (https://saves.mbi.ucla.edu/). Using these assessments, we firstly selected the best model from each algorithm and then compared them with each other. Since the qualities of the two best models were so close, they were further evaluated for folding stability via molecular dynamics (MD) simulations. Of note, before the quality assessments, the models subjected to MD simulations (see “further models” and “variant models” in Table 2) were processed with “RepairPDB” command of FoldX (v5.0) (Delgado et al., 2019; van Durme et al., 2011), which converts the saved coordinates to the PDB format of SAVES and applies a basic minimization procedure.

During the revision process of this article, a revolution in protein structure prediction was carried out in the Critical Assessment of protein Structure Prediction (CASP) (Pereira et al., 2021) by Google DeepMind. They developed a deep learning-based algorithm called AlphaFold2 (AF2) which makes de novo structure prediction with high accuracy (Jumper et al., 2021; Thornton, Laskowski & Borkakoti, 2021; Varadi et al., 2021). Because of its success, we compared our model with the AF2 model as an additional validation. AlphaFold Protein Structure Database (https://alphafold.ebi.ac.uk/entry/Q8TC57) (Varadi et al., 2021) was used to obtain AF2’s M1AP model.

We used PyMOL’s iterative superposition algorithm (Schrödinger & DeLano, 2020) for the structural alignment of the AF2 model with our model (see Fig. S1) and to compare the models generated by the same homology modeling algorithm (see Fig. S2).

Visual Molecular Dynamics (VMD v1.9.3) (Humphrey, Dalke & Schulten, 1996) was used to model the structures with the variants via the VMD Mutator Plugin.

System preparation

PDB2PQR (Dolinsky et al., 2004) was used to calculate the protonation states of amino acids at pH 7. Moreover, Glu43, Glu436 and Asp118 in the GalaxyWeb models (i.e., wild-type and all mutants) were protonated whereas His185, His238 and His375 were protonated in the Robetta-domain model. Models were solvated in a TIP3P water (Jorgensen et al., 1983) box. The systems were neutralized and ionized by adding 0.15 M KCl. The dimensions of the system were approximately $112.6\times 108.8\times 125.2$ ų in the xyz axes and the system contained ~146,000 atoms. VMD was used for the entire preparation procedure.

Simulation setup

NAMD 2.13-multicore-CUDA (Phillips et al., 2005) was used for MD simulations with the CHARMM36m all-atom force field (Huang et al., 2016). Simulations were performed under NPT ensemble with Langevin thermostat and Nosé-Hoover Langevin piston at 310 K and 1 atm (Feller et al., 1995; Davidchack, Handel & Tretyakov, 2009). The damping coefficient (gamma) of the Langevin thermostat was 1/ps. The oscillation period of the piston was around 100 fs whereas the oscillation decay time was 50 fs. A cutoff distance of 12 Å was used for the van der Waals interactions. The switching function starts at 10 Å and reaches zero at 14 Å. Long-range Coulomb interactions were calculated with particle-mash Ewald (Kolafa & Perram, 1992). A 2 fs integration time-step was used. The minimization and equilibration procedure was performed as follows: (1) 5,000 steps conjugate gradient algorithm and 1 ns equilibration of the solution by fixing the protein; (2) 5,000 steps of conjugate gradient minimization and 1 ns equilibration of the whole system by releasing all constraints except the ShakeH algorithm. This procedure provided sufficient relaxation for the system (Fig. S3). ShakeH of NAMD was applied during the entire simulation procedure to keep hydrogens constrained. Production simulations (Repeat 1) for each system, including wild-type and the variants, were performed for 500 ns. Production simulations were re-performed for 500 ns (Repeat 2) with different random number generator seeds to validate the reproducibility of the results and/or to obtain different possible conformations independent of the original simulations. In this way, the total production simulation time for each system was 1,000 ns. All simulation results are available at Zenodo (https://doi.org/10.5281/zenodo.5811977).

Structural and dynamical analyses

The “Generate” function of PDBsum (Laskowski, 2009) was used to predict the 2-D structure topology, clefts and tunnels of the model.

Alterations in a protein’s compactness, folding, solvent accessibility, stability, flexibility, hinges, local interactions (i.e., hydrogen bonds (H-bonds) and salt bridges within 10 Å of variant sites) and essential motions might influence functions, interactions and regulations of that protein. Therefore, we used in-house tcl scripts (available at https://github.com/ugerlevik/M1AP_analysis) and built-in VMD plugins to calculate backbone root-mean-square deviation (RMSD, related to folding), average C_α root-mean-square fluctuation (RMSF, indicates the flexibility), radius of gyration (Rg, inversely associated with compactness), solvent accessible surface area (SASA) and distances.

FoldX (v5.0) (Delgado et al., 2019) was used for stability analysis. ProDy (Bakan, Meireles & Bahar, 2011) was used for principal component analysis (PCA), which refers to conformational dynamics or essential motions related to the free energy landscape observed in the MD simulations. ProDy was also used for Gaussian network modeling to predict the hinge residues that play a key role in the essential motions of proteins.

Secondary structure analysis and related visualizations were performed by using the ggstride package (https://github.com/ugerlevik/ggstride/), which was developed for the first time to be used in this study. ggstride uses STRIDE (Frishman & Argos, 1995) for secondary structure assignment, bio3d (Grant et al., 2006) for structure and trajectory evaluations, and ggplot2 (Wickham, 2016) and ggpubr (Kassambara, 2020) for visualizations. This whole procedure was performed in R (v4.0.3) (R Development Core Team, 2020).

In PCA and secondary structure analysis, repeat 1 and 2 trajectories were merged while obtaining the observed conformations and secondary structure percentages. Only the last 75 ns of each trajectory was used in RMSF, PCA, secondary structure, and hinge analyzes as these parts had the most stable RMSD levels (Fig. 1) indicating that the structures had reached the energetic minimum.

Analysis results were visualized using the ggplot2, ggpubr and gg.gap (https://github.com/ChrisLou-bioinfo/gg.gap) packages in R, and the matplotlib package (Hunter, 2007) in python (v3.7.6) (van Rossum & Drake, 2009). Structural visualizations were performed by using VMD.

Evolutionary conservation

We observed that all variant sites except Ser50 were conserved among all orthologs examined (Fig. S4). Ser50 was also highly preserved by changing to a Cys only in the cat (Felis catus, UniProt ID: A0A337RXR8). This suggests that all these variants may have significant effects on the structure and function of meiosis 1-associated protein (M1AP). Therefore, all of them were further investigated using molecular dynamics (MD) simulations.

Molecular modeling and validation

Twenty-six models for the M1AP structure were built via the following algorithms with automatically selected templates: GalaxyWeb (five models), I-TASSER (five models), Phyre2 (one model), PRIMO (four models), RaptorX (one model), Robetta (five models) and Robetta-domain (five models). Although the properties regarding to the templates such as identity and coverage (Tables S1 and S2) indicate that the M1AP structure is difficult to model, it was possible to model the entire structure as the templates covered different parts of the M1AP and as utilizing the ab initio modeling. However, we required a comprehensive validation procedure as follows to proceed with the simulations reliably.

Firstly, we selected the best model from each algorithm by comparing their quality assessment scores (Table S3) to simplify benchmarking the models of the different algorithms. Remarkably, the models generated by the same algorithm were mostly similar to each other, except for the Robetta and Robetta domain models, which were more diverse than the models of the other tools (Fig. S2). As a result, model 1 from GalaxyWeb, I-TASSER, PRIMO and Robetta-domain, and model 5 from Robetta came forward for further benchmarking. Furthermore, among the single models of RaptorX and Phyre2 and the selected models of other tools, the top three models were from GalaxyWeb, Robetta and Robetta-domain (see “Initial Models” in Table 2). Since the GalaxyWeb model had the highest Ramachandran plot (PROCHECK) scores and the Robetta-domain model had slightly better scores than the Robetta model, we retained the GalaxyWeb and Robetta-domain models to continue with them as the two best models.

Instead of basic and static algorithms, we compared the two best models based on the observation of their folding stability over 100 ns of molecular dynamics (MD) simulations. We used the backbone root-mean-square deviation (RMSD) as an indicator of the stability of predicted protein folding. As shown in Fig. 2A, the GalaxyWeb model was quite stable after a short time whereas the Robetta-domain model continued to deviate throughout the entire simulation period. In addition, the quality scores of the models at the end of the 100 ns simulations differed only slightly from the initial structures, as expected (see “Further Models” in Table 2). Based on the RMSD results, we finally chose the GalaxyWeb model to examine the effects of variants. This M1AP model is available in the PMDB Protein Model DataBase (http://srv00.recas.ba.infn.it/PMDB/) (Castrignanò et al., 2006) under the accession identifier PM0083500.

Since the AlphaFold2 (AF2) revolutionized the field of de novo protein modeling during the revision process of this article (Jumper et al., 2021; Thornton, Laskowski & Borkakoti, 2021; Varadi et al., 2021), we compared our final model with the AF2 model as an additional validation. As shown in Fig. S1, our model and AF2 model were quite similar (RMSD = 8.42 Å). In addition, the AF2 model has similar level of quality (i.e., better VERIFY3D and PROVE but lower ERRAT, PROCHECK and WHATCHECK scores) with our final model (see “Comparison Model” in Table 2). Because of this independently obtained resemblance of our GalaxyWeb model to this revolutionary de novo structure, we might infer that it is one of the best structures that can be computationally obtained by present knowledge.

Variant models were evaluated to see if they retained quality. As shown in “Variant Models” in Table 2, their qualities were quite similar to those of the wild-type near the end of the simulations (i.e., 450^th ns), pointing out the validity of the simulation results.

Structural details

The final M1AP model has a globular shape (Fig. 2B) and involves 13 α-helices and 11 β-sheets. The detailed topology including the secondary structures generated by PDBsum is shown in Fig. 2C. Moreover, the model has a large cleft (volume: 11,448 ų) shown in Fig. 2D. Upon observations of the absence of M1AP in the nucleus and impaired meiosis in the M1AP deficiency, Arango et al. (2013) theorized that M1AP might indirectly regulate meiosis progression by binding and activating a cytoplasmic protein that would later translocate into the nucleus. If this suggestion is correct, this large cleft in M1AP would be the most probable binding site for this theorized protein. Other clefts are not shown because they were significantly smaller (i.e., 3272 ų, 2207 ų, and 1507 ų) and less likely to be suitable for ligand binding (Laskowski et al., 1996).

The locations of missense variants are shown on the tertiary structure of M1AP in Fig. 2D (also visible in the secondary structure topology of M1AP in Fig. 2C). Moreover, Ser50 and Leu430 are in α-helices; Gly317 is in a loop close to an α-helix; Pro389 is at the edge of a β-sheet; and Arg266 is in a surface-exposed loop. On the other hand, Pro389 and Leu430 are part of the largest cleft while Gly317 is close to it (Fig. 2D).

Impacts of the variants

Folding, compactness and solvent accessibility

The most stable part of the RMSD was the last 75 ns of each 500 ns trajectory, indicating that the simulations reached equilibrium (Fig. 1). Taking this into account, backbone deviation was lower in variants p.Ser50Pro (both repeats) and p.Arg266Gln (only repeat 1) and slightly higher in p.Gly317Arg (only repeat 1) than in wild-type whereas there were no remarkable changes in other variants (Fig. 1). Of note, a rapid and sharp increase was observed in the p.Leu430Pro variant before gradually decreasing and stabilizing. Combined with the radius of gyration (Rg) pattern (Fig. 3), it can be inferred that the structure partially unfolded and refolded during this time.

Rg was slightly lower in variants p.Arg266Gln (both repeats), p.Gly317Arg (only repeat 2) and p.Leu430Pro (both repeats), indicating that the protein was more compact upon these amino acid changes (Fig. 3). Rg of other variants did not differ from the wild-type. In any variant, we did not observe any critical solvent accessible solvent area (SASA) differences with wild-type (Fig. 4), which might be expected upon compactness changes.

Flexibility, hinge and essential dynamics

The conformational spaces observed in all variants were significantly different from the wild-type (Fig. 5), indicating altered essential protein dynamics and energetic minima of backbone conformations. In line with this observation, there were critical root-mean-square fluctuation (RMSF, indicating residue flexibility) changes in at least one repeat of each variant as pointed out by the yellow stars in Fig. 6, and hinge site (zero flexibility regions) differences as indicated by yellow and red colors in Table S4.

Secondary structures

Ser50 and Leu430 were found in α-helices (Figs. 2C and 2D) and changed to Pro, which is one of the amino acids known as secondary structure breakers. Consistent with this expectation, both variants disrupted their α-helices, making them turn and 3(10)-helix structures (Fig. 7), which are less stable than α-helices. In p.Arg266Gln, the variant region was more stable as α-helix and turn rather than a coil observed in wild-type (Fig. 7). Furthermore, these were the most outstanding secondary structure changes at variant sites. However, there were less significant changes in other parts of the protein as well (Fig. S9).