# Deciphering the Binding of Salicylic Acid to Arabidopsis thaliana Chloroplastic GAPDH-A1

Igor Pokotylo, Denis Hellal, Tahar Bouceba, Miguel Hernandez-Martinez, Volodymyr Kravets, Luis Leitao, Christophe Espinasse, Isabelle Kleiner, Eric Ruelland

► To cite this version:

Igor Pokotylo, Denis Hellal, Tahar Bouceba, Miguel Hernandez-Martinez, Volodymyr Kravets, et al.. Deciphering the Binding of Salicylic Acid to Arabidopsis thaliana Chloroplastic GAPDH-A1. International Journal of Molecular Sciences, 2020, 21 (13), pp.4678. 10.3390/ijms21134678. hal-03017984

HAL Id: hal-03017984

https://hal.science/hal-03017984

Submitted on 23 Nov 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

Deciphering the Binding of Salicylic Acid to
*Arabidopsis thaliana* Chloroplastic GAPDH-A1

Igor Pokotylo <fnref n="1" />,<fnref n="2" />, Denis Hellal <fnref n="1" />, Tahar Bouceba <fnref n="3" />, Miguel Hernandez-Martinez <fnref n="1" />>D,
Volodymyr Kravets <fnref n="2" />, Luis Leitao <fnref n="1" />, Christophe Espinasse <fnref n="1" />, Isabelle Kleiner <fnref n="4" /> and
Eric Ruelland <fnref n="1" />,*D

<figure><img src="image_3.png" /></figure>

check for
updates

Received: 27 May 2020; Accepted: 28 June 2020; Published: 30 June 2020

**Abstract:** Salicylic acid (SA) has an essential role in the responses of plants to pathogens. SA initiates defence signalling via binding to proteins. NPR1 is a transcriptional co-activator and a key target of SA binding. Many other proteins have recently been shown to bind SA. Amongst these proteins are important enzymes of primary metabolism. This fact could stand behind SA's ability to control energy fluxes in stressed plants. Nevertheless, only sparse information exists on the role and mechanisms of such binding. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was previously demonstrated to bind SA both in human and plants. Here, we detail that the A1 isomer of chloroplastic glyceraldehyde 3-phosphate dehydrogenase (GAPA1) from *Arabidopsis thaliana* binds SA with a *K*<sub>D</sub> of 16.7 nM, as shown in surface plasmon resonance experiments. Besides, we show that SA inhibits its GAPDH activity in *vitro*. To gain some insight into the underlying molecular interactions and binding mechanism, we combined in silico molecular docking experiments and molecular dynamics simulations on the free protein and protein-ligand complex. The molecular docking analysis yielded to the identification of two putative binding pockets for SA. A simulation in water of the complex between SA and the protein allowed us to determine that only one pocket—a surface cavity around Asn35—would efficiently bind SA in the presence of solvent. In silico mutagenesis and simulations of the ligand/protein complexes pointed to the importance of Asn35 and Arg81 in the binding of SA to GAPA1. The importance of this is further supported through experimental biochemical assays. Indeed, mutating GAPA1 Asn35 into Gly or Arg81 into Leu strongly diminished the ability of the enzyme to bind SA. The very same cavity is responsible for the NADP+ binding to GAPA1. More precisely, modelling suggests that SA binds to the very site where the pyrimidine group of the cofactor fits. NADH inhibited in a dose-response manner the binding of SA to GAPA1, validating our data.

**Keywords:** salicylic acid; glyceraldehyde 3-phosphate dehydrogenase; molecular dynamics; molecular docking; protein ligand interaction; surface plasmon resonance; biacore

<footnote n="1">¹ IEES-Paris (UMR 7618)—Institut D’écologie et des Sciences de L’environnement de Paris; CNRS UMR 7583, Université Paris-Est Créteil, Sorbonne Université, F-94010 Paris, France; pokotylo@bpci.kiev.ua (I.P.); denishellal@gmail.com (D.H.); miangher2010@gmail.com (M.H.-M.); luis.leitao@u-pec.fr (L.L.); christophe.espinasse@u-pec.fr (C.E.)</footnote>

<footnote n="2">² V.P. Kukhar Institute of Bioorganic Chemistry and Petrochemistry, National Academy of Sciences of Ukraine, 02094 Kyiv, Ukraine; kravets@bpci.kiev.ua</footnote>

<footnote n="3">³ Plateforme D'interactions Moléculaires, CNRS-FR3631; Institut de Biologie Paris Seine (IBPS), Sorbonne Université, Cedex 05, F-75252 Paris, France; tahar.bouceba@upmc.fr</footnote>

<footnote n="4">⁴ LISA (UMR 7583)—Laboratoire Interuniversitaire des Systèmes Atmosphériques (LISA), CNRS UMR 7583, Université Paris-Est Créteil, Université de Paris, Institut Pierre Simon Laplace (IPSL), 61 Avenue du Générale de Gaulle, F-94010 Créteil, France; kleiner@lisa.u-pec.fr</footnote>

<footnote>* Correspondence: eric.ruelland@u-pec.fr</footnote>
# 1. Introduction

Salicylic acid (SA, 2-hydroxybenzoic acid) is an important phytohormone that controls defence reactions in plants challenged by pathogens [1,2]. SA is produced in the cells located at the site of the pathogen attack. In these cells, SA initiates downstream responses aiming at counteracting the attack. Besides, the methylation of SA into methyl salicylate (MeSA) by a SA methyl transferase increases its volatility and cellular permeability. MeSA could serve as a chemical signal travelling from the site of the infection to distant tissues. In these distant sites, MeSA is de-methylated back to SA, which then contributes to establishing a systemic immunity (systemic acquired resistance, SAR) [3,4]. SA is also a part of some abiotic stress reactions in plants, e.g., drought responses in oat [5] or salinity tolerance in *Arabidopsis thaliana* [6].

SA also has a role in animals, where salicylate-derived compounds are indeed being extensively used in medicine. Aspirin is acetylsalicylate. When administered orally, aspirin is rapidly converted into SA, which has a much longer half life as compared to its acetylated counterpart [7].

Signalling pathways initiated by SA in the cells are starting to be understood [8]. However, while SA has an essential role, either in the physiology of plants or both in human and animal health, only sparse information exists on the role of proteins that bind SA and on the mechanisms of these SA-protein interactions [2].

Interestingly, it seems that SA recognition and signalling does not implicate membrane receptors that would trigger cascades of protein kinases, leading to activation or the inhibition of target enzymes. This would represent a classical hormone pathway described in mammals and for plant brassinosteroids. On the contrary, it appears that SA directly binds to cytosolic proteins, potentially modulating their functions or enzymatic activity. SA binds to NPR1, which is an essential transcription co-regulator of the SA pathway. SA has a role in NPR1 monomerisation [9]. NPR1 monomers shuttle from the cytosol to the nucleus where they interact with TGA transcription factors, thus allowing the induction of SA-responsive genes in plants [8]. NPR1 is a critical protein component of the SA pathway and a great deal of SA sensitivity is lost in *npr1* plants [10]. Nevertheless, the binding of SA to NPR1 cannot explain the whole range of SA responses. First, some responses to SA seem to be NPR1-independent [11–13]. Recently, SA-induced attenuation of root growth in Arabidopsis was attributed to SA binding to A subunits of protein phosphatase 2 (PP2A) and the inhibition of its activity acting on PIN auxin efflux carriers. Such an effect of SA is well preserved in *npr1* plants [14]. Second, even within the NPR1-dependent pathway, SA binding to NPR1 does not explain it all. The monomerisation of NPR1 occurs through the reduction of disulphide bridges in the NPR1 homotetramer; binding to SA alone is not enough to monomerise NPR1 [9]. Moreover, NPR1 monomerisation could altogether be SA-independent. For instance, endoplasmic reticulum (ER) stress promoted NPR1 nuclear translocation, which could be only due to its monomerisation. During ER stress, SA is not accumulated but a redox potential of the cytosol is reduced [15]; this is the same effect observed following SA treatment [16]. This suggests that the role of SA in NPR1 monomerisation could be due to targeting protein regulators of redox status and not due to NPR1 binding *per se*. The cytosolic thioredoxins TRX-h3 and TRX-h5 were pulled down with the N-terminal part of NPR1 and were suggested to act in NPR1 monomerisation [17]. The regulation of their activity by SA is not confirmed. Finally, we have shown that SA rapidly activates phosphatidylinositol-4-kinases [18] and phospholipase D [19]. These activations are necessary for the triggering of the NPR1-dependent pathway [18,19]. Here again, SA response needs to occur through binding proteins that will trigger these lipid-signalling pathways upstream of NPR1.

Proteins other than NPR1 that bind SA have been found. The first described salicylic acid binding protein (SABP1) was a catalase, which was identified in tobacco [20]. SABP2 is a methysalicylate esterase [21,22], and SABP3 is a beta carbonic anhydrase [23]. With the exception of SABP2, their exact roles in the SA signalling pathway were not documented. SA has also been shown to bind NPR1 paralogs [24,25]. Recently, by high-throughput biochemical screenings, nearly 100 soluble proteins were proposed to be able to bind SA in *Arabidopsis thaliana* plant cells [26]. Such a large number of proteins that potentially could bind SA is surprising. Indeed, for many hormones, only a limited
number of proteins playing the role of receptors are usually found. Moreover, in this list of proteins, essential enzymes for the primary metabolism exist such as the small unit of the RUBISCO protein (AtRBCS1A, the key enzyme of the CO₂ fixation during photosynthesis) and several isomers of *Arabidopsis thaliana* chloroplastic GAPDH that are important enzymes of photosynthesis: GAPA1 and GAPA2 [27]. Therefore, the understanding of the mechanism of the interaction between a protein and SA is more than ever a relevant biological question.

In humans, the mechanisms of salicylate or aspirin actions have also proven to be much more complex than initially suggested (having a role in cyclooxygenase inhibition). A number of protein targets of aspirin were identified [28]. SA also could bind human cytosolic GAPDH, acting in glycolysis [29]. Even though SA is not a hormone in animals, evolutionarily, this can be explained by SA uptake by animals from plant food sources.

Despite our knowledge of proteins that bind SA, we still know little about the mechanisms of such binding. For instance, in plants, it was shown that SA binds and inhibits catalase SABP1 [30] or thimet oligopeptidases TOP1 and TOP2 in vitro [31]. How this is happening on the molecular level is still unclear.

To elucidate the binding motifs and mechanism at an atomistic scale, combined methods of in silico simulations (molecular dynamics) with in vitro binding assays were applied. We started this study with a glyceraldehyde 3-phosphate dehydrogenase (GAPA1), which is a member of the evolutionary conserved glyceraldehyde 3-phosphate dehydrogenase (GAPDH) family. Our approach allowed us to identify a surface cavity, structured by Asn35, as the one allowing the binding of SA. Asn35 and Arg81 were shown to be important for such binding. Interestingly, this site is the very same site where the adenine part of the cofactor NAD(P) is accepted by GAPA1.

The approach we used can be applicable to other SA target proteins, either of plant or animal origin.

## 2. Results

### 2.1. GAPA1 Binds SA with High Affinity

GAPA1 was produced in *E. coli* and purified to more than 90% (Figure 1). Then, we studied its interaction with SA via surface plasmon resonance (SPR). In this technique, an aminoethyl derivative of SA—3-aminoethyl salicylic acid (3-AESA)—is immobilised on the sensor surface and serves as bait. A running buffer containing GAPA1 is flown over 3-AESA to detect binding in real time. GAPA1 bound 3-AESA in a concentration-dependent manner. Representative sensorgrams are available in Supplementary Figure S1A. These binding experiments were fitted to a 1:1 Langmuir model. From this model, an equilibrium dissociation constant ($K_D$) could be calculated. The GAPA1-3-AESA interaction yielded to a mean $K_D$ of 16.7 ± 5.7 nM in three independent assays. From the Langmuir modelling, responses at equilibrium (Req) could be extracted and plotted against the concentration of protein (Figure 2A). Note that no binding was observed when bovine serum albumin (BSA; used as a negative control) was flown over 3-AESA (Supplementary Figure S2).

Following the observation of a GAPA1 binding to 3-AESA, we then tested if this interaction was relevant in regard to the interaction with native SA. To do so, we have performed competition tests in the presence of an excess of SA. Relative responses were extracted just after the beginning of the dissociation phase of the sensorgrams (Supplementary Figure S1B,C). GAPA1 association with immobilised 3-AESA was inhibited by the presence of SA in a running buffer in a concentration-dependent manner (Figure 2B). This demonstrated that SA competes with 3-AESA for GAPA1 binding (i.e., binding occurs at the same site/sites). The presence of 4-hydroxybenzoic acid (4-HBA), a biologically inactive structural isomer of SA, had led to a much weaker competition (Figure 2B).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Production and purification of recombinant glyceraldehyde 3-phosphate dehydrogenase (GAPA1). In (A), a detection of GAPA1 with anti-6xHis antibodies in a soluble fraction of *E. coli* lysate via the peroxidase activity of secondary antibodies. In (B), a SDS-PAGE of GAPA1 purified from a soluble fraction of *E. coli* lysate using Ni-NTA resin. The protein band was visualised by stain-free imaging (Bio-Rad). Shading at the bottom of the track is due to the loading buffer pigments.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Assessment of GAPA1 binding to immobilised 3-aminoethyl salicylic acid (3-AESA) in surface plasmon resonance (SPR) assay. In (A), a range of GAPA1 concentrations were tested for their ability to bind immobilised 3-AESA. K_D of this interaction was calculated by fitting the binding profiles to a 1:1 Langmuir interaction model (Supplementary Figure S1A). Binding experiments were performed 3 times and led to a mean K_D value of 16.7 ± 5 nM. The Langmuir modelling allowed extracting the response at equilibrium (Req) values for each tested protein concentration. Req values were plotted against GAPA1 concentrations. In (B), GAPA1 ability to bind 3-AESA was assessed in the presence of an excess of SA or 4-HBA. For that, the concentration of GAPA1 was set at 200 nM. Samples were pre-incubated with an indicated concentration of SA on ice for 1 h prior to the binding test. Trend lines are generated for visualisation purposes.</figcaption></figure>
## 2.2. SA Inhibits GAPA1 Catalytic Activity in Vitro

Then, we wanted to analyse the effect of SA on GAPDH enzymatic activity. GAPA1 was pre-incubated with SA on ice for 1 h prior to in vitro enzymatic assays. SA inhibited GAPA1 activity in a concentration-dependent manner (Figure 3A). Incubation with an equivalent concentration of 4-HBA did not lead to an inhibitory effect: at 2mM, SA produced 28% inhibition of GAPDH activity of GAPA1, while 2 mM of 4-HBA produced an inhibition of only 4% (statistically insignificant) (Figure 3B).

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) activity (NAD-dependent) of GAPA1 pre-incubated with SA or 4-hydroxybenzoic acid (4-HBA) on ice for 1 h. Graphs display averages of 3 replicates. Error bars show SD. Statistical differences were confirmed by ANOVA (panel A) and by a two-tailed Student's t test (panel B): \* (α = 0.05), \*\*\* (α = 0.001).</figcaption></figure>

## 2.3. In Silico Modelling of SA Interaction with GAPA1

Our data strongly indicate that SA binds GAPA1, but these biological assays do not allow for a detailed understanding of the ongoing phenomena at the molecular level. Here, in silico techniques such as molecular dynamics simulations are the methods of choice to elucidate the underlying interaction mechanism at an atomistic scale.

We extracted the structure of one GAPA1 monomer from the structure of the wild-type protein AtGAPDHA1 crystallised as a tetramer (PDB 3K2B; Supplementary Figure S3). Note that even though NADP is the preferred cofactor of the chloroplastic GAPA1, the authors crystallised it with NAD [32].

This initial structure was relaxed in water for 600 nsec using GROMACS package v5.0.7. The root mean square deviation (RMSD) of atomic positions, which is a measure of global backbone deviation, was calculated using the structure of the protein before the simulation started (i.e., after the pressure equilibration step) as the reference. The RMSD presents a quick increase that attains a plateau after 40 nsec; yet, a big decrease of about 0.2 nm is observed at 200 nsec of simulation; then, a sharp increase in RMSD occurs at around 410 nsec (Figure 4). A root mean square fluctuation (RMSF) was calculated between 0 and 75 nsec: the rapid change of the RMSF occurs in the vicinity of the residue Ala200 that corresponds to an unstructured loop (Supplementary Figure S4A). RMSF was also calculated for 190–250 nsec and 400–500 nsec simulation spans: the decrease in RMSD at 200 nsec and the subsequent increase are also due to fluctuation around Ala200 (Supplementary Figure S4A).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Relaxation of the protein structure in water. Root mean square deviations (RMSD) are shown for GAPA1 (GROMACS package v5.0.7). In RMSD, each structure from a trajectory is compared to a reference structure; here, the reference used is the structure of the protein after the pressure (NPT) equilibration step. The time frames corresponding to the 5 clusters are shown in colour. Clustering was performed for the simulation time 40 nsec to 600 nsec.</figcaption></figure>

Then, we performed a clustering of the structures obtained during the relaxation from 40 to 600 nsec with a sampling of 390 psec. We obtained 5 clusters (Figure 4). Clusters 1 to 3 comprise 73% of the simulation structures. Cluster 1 mostly corresponds to structures of the first plateau (40–190 nsec); cluster 3 mostly corresponds to structures of the 450–600 nsec plateau; cluster 2 is comprised of structures within the valley between 200 and 400 nsec of the simulation. For each cluster, an average structure was generated. By overlapping the average structures for clusters 1 to 3, we could visualise that the main difference between structures is due to the loop comprising Ala200 (Supplementary Figure S4B). Note that the second increase in RMSD leading to the second plateau (450–600 nsec) does not indicate reverting to structures similar to that of the first plateau (40–190 nsec). Indeed, when RMSD was calculated using the average structure 1 as the reference (Supplementary Figure S4C), we could see that structures in the second plateau (450–600 nsec) are in fact more distant to the structures of the first plateau (that are obviously the closest to average structure 1) than the structures of the valley (200–400 nsec).

### 2.3.1. Preparation of the Ligand: Calculation of Atomic Charge of the Ligands

Molecular geometries of SA have been optimised using the *ab initio* method at MP2/6-311++G** level of theory using GAUSSIAN09 (https://gaussian.com). The SA PDB file used in our docking can be found as Supplementary File 1.

### 2.3.2. Determination of SA-Binding Pockets

AUTODOCK VINA [33,34] was used to perform the docking of SA in the average structures of clusters 1 to 3. Docking was performed without water molecules and with flexible side chains. As we did not know the binding pocket, a blind docking was performed enclosing the whole protein into the grid box. Standard default parameter settings were used. When considering these 10 best docking scores in the average structures of clusters 1 to 3, the SA was located in pockets labelled A to H (Table 1). Figure 5 illustrates those 10 best docking scores of SA in the average structure of cluster 3.
<table><caption>**Table 1.** Docking of salicylic acid (SA) in the average structures of clusters 1, 2 and 3. Docking was performed with AUTODOCK VINA [33,34] implemented in PyMOL. Spacing was set at 0.375 for grid definition.</caption><thead><tr><th colspan="3">Average Structure Cluster 1</th><th colspan="3">Average Structure Cluster 2</th><th colspan="3">Average Structure Cluster 3</th></tr><tr><th>Rank <sup>a</sup></th><th>Score <sup>b</sup></th><th>Pocket <sup>c</sup></th><th>Rank <sup>a</sup></th><th>Score <sup>b</sup></th><th>Pocket <sup>c</sup></th><th>Rank <sup>a</sup></th><th>Score <sup>b</sup></th><th>Pocket <sup>c</sup></th></tr></thead><tbody><tr><td>1</td><td>-5.9</td><td>A</td><td>1</td><td>-6</td><td>F</td><td>1</td><td>-6.3</td><td>F</td></tr><tr><td>2</td><td>-5.6</td><td>A</td><td>2</td><td>-5.9</td><td>F</td><td>2</td><td>-6.2</td><td>F</td></tr><tr><td>3</td><td>-5.5</td><td>A</td><td>3</td><td>-5.8</td><td>F</td><td>3</td><td>-5.7</td><td>G</td></tr><tr><td>4</td><td>-5.4</td><td>A</td><td>4</td><td>-5.7</td><td>F</td><td>4</td><td>-5.5</td><td>F</td></tr><tr><td>5</td><td>-5.3</td><td>B</td><td>5</td><td>-5.5</td><td>F</td><td>5</td><td>-5.5</td><td>F</td></tr><tr><td>6</td><td>-5.2</td><td>A</td><td>6</td><td>-5.4</td><td>F</td><td>6</td><td>-5.4</td><td>A</td></tr><tr><td>7</td><td>-5.2</td><td>C</td><td>7</td><td>-5.3</td><td>A</td><td>7</td><td>-5.4</td><td>F</td></tr><tr><td>8</td><td>-5.2</td><td>B</td><td>8</td><td>-5.3</td><td>F</td><td>8</td><td>-5.3</td><td>A</td></tr><tr><td>9</td><td>-4.9</td><td>D</td><td>9</td><td>-5.2</td><td>C</td><td>9</td><td>-5.3</td><td>A</td></tr><tr><td>10</td><td>-4.8</td><td>E</td><td>10</td><td>-5.2</td><td>A</td><td>10</td><td>-5.2</td><td>H</td></tr></tbody></table>

For the different protein structures, the docking poses are sorted from 1 to 10 (a) according to the calculated scores by AUTODOC VINA as a sum of intermolecular and intramolecular contributions (b). Each docking pose was visualised in the protein structure, thus allowing the naming of the pockets (c). Some docking poses are indeed in the same pockets.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Docking positions of SA in GAPA1. The 10 best scores docking poses are shown in the average structure of cluster 3. SA is in green and the protein is in magenta. Docking was calculated by AUTODOCK VINA [33,34] implemented in PyMOL. Structures were visualised using CHIMERA [35].</figcaption></figure>

From Table 1, even though the scores are sometimes close to each other, it appears that two cavities, A and F, combine the highest scores and the main occurrence from the 3 clusters. For the average structure of cluster 1, the most represented cavity is cavity A, for which the highest scores are obtained. For the average structure of clusters 2 and 3, cavity F is the most represented and with the highest scores. In these two average structures, cavity A is nevertheless present amongst the 10 best scores, at positions 7 and 10 for average cluster 2 and positions 6, 8, and 9 for average cluster 3.

The GAPA1 monomer structure can be broadly defined as two sub-globuli with each globule being structured by a beta sheet (Supplementary Figure S3). Pocket A is within the amino terminal globuli; it is situated at the top of the beta sheet; it involves residues of the beta sheet such as Asn35 and Asn8, and residues in the loops between the end of the beta sheets and helices, such as Arg81 and
Thr37. The upper side of cavity A is indeed capped by Arg81 and Thr37, while the bottom part of the cavity is comprised of Val102, Phe103, and Thr100. The inner part of the cavity is made of Asn35, Asp36, and Asn8 (Supplementary Figure S5A). Cavity F is at the interface of the two subparts of the peptide defined previously. It consists of Thr52 and Asp50 at one side and Ser290, Cys289, and Trp319 at the other side (Supplementary Figure S5B). When analysing the different docking possibilities of SA in cavity A (docking 1 to 4 and 6 in average structure of cluster 1), it appears that SA in cavity A is orientated in such a way that the negatively charged parts are towards the inside of the cavity, facing Asn35. Besides, the benzenic cycle is not necessary orientated the same way in all of the dockings. On the other hand, in pocket F (docking 1 to 6 and 8 in average structure of cluster 2), the benzenic cycle is always oriented parallel to the aromatic cycle of Trp319 (Supplementary Figure S6)). This suggests the importance of the pi-pi interactions in SA binding to cavity F.

### 2.3.3. Molecular Dynamics Simulations of the Complexes between SA and the Native GAPA1

For both cavities A and F, we performed five 20 ns simulations of the protein/SA complex in water with GROMACS package v5.0.7 to make sure that the complex is stable. We calculated the distance of SA to Asn80 and to Trp319 for pocket A and pocket F, respectively. We considered that SA is no longer bound to the pocket if the distance to the reference residue exceeded 1 nm in the course of the 20 nsec simulation. Interestingly, out of 5 simulations, SA drifted away from pocket F in 4 of these simulations. On the contrary, SA remained in the pocket A for this period of time in 4 out of 5 20 nsec simulations. In these simulations, the minimal distance from SA to Asn80 fluctuated in the range of 0.23 nm to 0.91 nm. Representative examples of stay/leave behaviour of SA, as shown on distance graphs, can be seen at Figure 6A,B. In pocket A, at the beginning of the simulation, SA was positioned with the -OH and -COOH groups towards the inside of the pocket (Figure 6C), pointing towards Asn35, thus producing H bonds. Throughout the simulation, SA made a full spin so that these groups are positioned towards the exit of the pocket, interacting with Arg81 (Figure 6D). H-bonds and hydrophobic interactions were calculated using LigPlot⁺ v2.2 software. Throughout the 10 last nanoseconds of the complex simulation, the four residues calculated to be the closest to SA are Asn8, Asn35, Arg81, and Phe103 (Figure 6E).

### 2.3.4. Molecular Dynamics Simulations of the Complexes between SA and Mutated AtGAPA1

In silico mutagenesis were performed on selected residues of pocket A. Similarly to native proteins (see above), 5 simulations were performed for each mutation. SA was considered to leave the cavity if the distance to Asn80 was higher than 1 nm (Figure 7).

Arginine is a residue with a positively charged lateral chain. When Arg81 was mutated into Ala or Leu, i.e., residues with hydrophobic side chains, SA left the cavity in 4 out of the 5 simulations in either case. When Arg81 was mutated into Gln, a residue with a polar but uncharged lateral chain, SA left the cavity in 3 out of the 5 simulations. It is also the case when Arg81 is mutated into Gly, for which the lateral chain is a H. This shows the importance of the positive charged lateral chain of Arg81 in the interaction with SA. When Arg81 was mutated into Lys, another positively charged residue, SA left the cavity in only 2 out of 5 simulations. We could observe that in the Arg81Lys mutant, the SA position in the cavity is also stabilised by producing an H-bond with Thr37 (Supplementary Figure S7A). This is also the case for some simulations in Arg81Ala and Arg81Gly (Supplementary Figure S7B,C). Therefore, Thr37, a residue that is not in close contact with SA in native pocket A (Figure 7E), could have some role in the docking of SA, especially in the case of Arg81 substitution. In the Thr37Ser and Thr37Cys mutants, SA left the pocket in 3 out of 5 simulations. Yet, when Thr37 is mutated into Gly for which the lateral chain is H, SA always stayed in the pocket. In this mutant, SA produced an H-bond with Arg81, similarly to the wild-type (WT) protein (Supplementary Figure S7D). Clearly, the role of Thr37 in the docking of SA is lesser to that of Arg81.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Molecular dynamic simulations of the complex between GAPA1 and SA. In (A), the distance between SA and Asn80 was calculated throughout a course of 20 nsec simulation; the starting docking corresponds to the best scored docking of SA in the average structure of cluster 1 (SA in pocket A). In (B), the distance between SA and Trp319 was calculated throughout a course of 20 nsec simulation; the starting docking corresponds to the best scored docking of SA in the average structure of cluster 3 (SA in pocket F). In (C), a position of SA is shown immediately after the temperature and pressure equilibrations using GROMACS package v5.0.7. In (D), a position of SA in the pocket is shown after a 20 nsec simulation in water. The starting docking used to generate panels (C) and (D) corresponds to the best scored docking of SA to the pocket A in the average structure of cluster 1. Structures were visualised using CHIMERA [35] and LigPlot+ [36]. H-bonds and distances (in nm) are shown as dashed green lines. Red spoked arcs represent protein residues making nonbonded contacts with the ligand. In (E), the distances between SA and each residue of GAPA1 were calculated during the 10 last nanoseconds of the simulation of SA docked to pocket A of the average structure of cluster 1.</figcaption></figure>
Asn35 was shown to be one of the closest residues to SA (Figure 6E). Asparagine has a polar yet uncharged lateral chain. Asn35 was mutated into Ala and Leu, which are residues with apolar lateral chains, or into Asp, a polar residue with a negatively charged chain, or Gly, for which the lateral chain is H. In all cases, this led to SA leaving the pocket in at least 4 out of the 5 simulations. As explained above, during the course of simulation in the WT protein, SA goes from a position where its -OH and -COOH groups are directed towards Asn35 (towards the inside of the pocket) to a position where they are positioned towards Arg81, i.e., to the outside of the pocket (Figure 6C,D). Hence, the importance of Asn35 is not to directly interact with SA, but more likely to structure the pocket. Indeed, Asn35 is shown to establish H-bonds with the backbone chain between Arg81 and Asn80 and Asn80 and Ser79; Asn35 is also calculated to establish a bond with Asn8 (data not shown). Asn35 is clearly at the interface of different peptide sections that constitute the pocket.

When Asn8 was mutated into Leu, SA left the pocket in 3 out of 5 simulations. Here again, SA could be stabilised by producing an H-bond with Thr37 (Supplementary Figure S7E). Phe103 defined one part of the entrance to the pocket. When Phe103 was mutated into Leu, SA left the pocket in 4 out of 5 simulations. This suggests that pi-pi interactions between the aromatic rings of Phe and SA could be of importance for retaining the latter in the cavity.

The residue 36 of GAPA1 is an aspartate, which is a polar negatively charged residue. Since it is near Asn35, we tested its possible role. When Asp36 was mutated into Ala or Gly, SA stayed in the pocket 4 and 3 times out of 5, respectively (Figure 7).

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Twenty nsec simulations of the complexes between SA and mutated GAPA1. SA was docked in pocket A before simulation started. This docking of SA corresponds to the docking scored 1 in Table 1 for the average structure of cluster 1. The in silico mutations are performed on the average structure of cluster 1. Simulations were run 5 times from 5 independent energy minimisation steps. Criteria for assuming ligand leaving the pocket: distance to Asn80 greater than 1 nm.</figcaption></figure>

### 2.3.5. SA is Docked to the NADPH Binding Pocket

We wanted to know whether the residues we identified are conserved amongst GAPDHs. We retrieved GAPDH sequences from animals and plants; for plants ones, we considered NAD-dependent ones and NADP-dependent ones as well (Figure 8). We could see that Asn35 and Arg81 are conserved in the assessed proteins. To what functions do they correspond? To understand
that, we overlapped the crystal structure of GAPDH-A1 (PDB 3K2B) with the structure obtained at 20 nsec of the SA/WT protein complex simulation. In their structural work [32], the authors crystallised GAPDHA1 in the presence of NAD even though NADP is the preferred cofactor of this enzyme. Interestingly, SA overlaps with the adenine group of the NAD cofactor (Figure 9).

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Alignment of different GAPDHs. Regions around Asn35 and Arg81 (marked with red triangles) of AtGAPA1 are shown. Plant NADP-dependent chloroplastic GAPDHs that are involved in photosynthesis are highlighted in green. AtGAPCP2 and AtGAPCP1 are NAD-dependent enzymes that function in chloroplasts. Animal and other plant GAPDHs are the cytosolic ones, while NAD-dependent ones are involved in glycolysis. Alignment was performed by ClustalX2. Conservation score is depicted using standard colour scheme of ESPrict [37]. Residues are shown on red background if strictly conserved. At, *Arabidopsis thaliana*; Nt, *Nicotiana tabacum*; Zm, *Zea mays*; So, *Spinacia oleracea*; Os, *Oryza sativa*; Pp, *Physcomitrella patens*; Ph, *Petunia hybrid*; Hs, *Homo sapiens*; Mm, *Mus musculus*; Oc, *Oryctolagus cuniculus*; Dm, *Drosophila melanogaster*.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** The SA binding pocket is the pocket where the adenine part of the NAD(P) cofactor fits. Structure in magenta: the protein structure in complex with SA after a 20 nsec simulation (such as in Figure 6A); structure in blue: the crystallographic structure corresponding to PDB 3K2B; orange: NAD; green: SA. (A) Overlapping between 3K2B and our simulation. (B) Close-up of A. Structures were visualised using CHIMERA [35].</figcaption></figure>
## 2.4. Testing the Ability of Point-Mutated GAPA1 to Bind SA

Based on in silico data, a cavity around Asn35 and involving Arg81 (cavity A) was suggested to be important for SA binding in GAPA1. Pre-selected residues of this cavity (Arg81 and Asn35) were validated by biochemical assays. Point-mutated GAPA1 proteins were produced and compared for their ability to bind SA in a similar SPR technique used for WT GAPA1. This allowed us to test if in vitro assays are in agreement with our theoretical predictions.

In GAPA1 Arg81Leu, the ability to bind 3-AESA (as shown by relative responses) was diminished by 55%. Whereas in GAPA Asn35Gly, it was diminished by 66% (Figure 10), suggesting that these amino acid residues contribute to SA binding to GAPA1. No complete loss of ability to bind 3-AESA was observed in these mutated proteins. This could be due to the fact that (unlike simulations obtained in silico), singular point mutation does not completely disrupt interaction with SA. Due to the problems with protein stability, we could not produce an Arg81Leu-Asn35Gly double-mutated GAPA1.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison of the abilities of wild-type (WT) GAPA1 (red line) and point-mutated GAPA1 Asn35Gly (pink line) and GAPA1 Arg81Leu (green line) to bind immobilised 3-AESA in SPR assays. All proteins were injected in 50 nM concentration. In these sensorgrams, the signal from the mock-coupled surface was subtracted. The response values of report points (*) after the beginning of the dissociation phase were extracted. The corresponding relative responses are indicated in Resonance Units (RU). The report points are indicated by *.</figcaption></figure>

In the docking calculations, we had also identified a second potential SA binding site (cavity F) implicating Trp319. This cavity had a much lower score of retaining SA in molecular dynamic simulations. Indeed, we found out that the ability of GAPA1 Trp319Leu to bind 3-AESA was not affected (Supplementary Figure S8).

## 2.5. *Ability of GAPA1 to Bind SA is Inhibited by NADH*

We have obtained preliminary data indicating that SA binding to GAPA1 could overlap with NAD(P) cofactor (Figure 9). In order to confirm such an observation, we have tested the ability of GAPA1 to bind immobilised 3-AESA in the presence of NADH. We used NADH as a competing molecule, since GAPA1 was crystallised with it [32]. Indeed, we have found that the pre-incubation of GAPA1 on ice with NADH reduced the level of GAPA1 association with immobilised 3-AESA (Figure 11). The competition was observed as soon as 2 µM NADH.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** The ability of WT GAPA1 to bind immobilised 3-AESA is inhibited in the presence of NADH. GAPA1 concentration was set at 200 nM concentration. The response values of report points (*) after the beginning of the dissociation phase were extracted. The corresponding relative responses are indicated in Resonance Units (RU). The report points are indicated by *. In these sensorgrams, the signal from the mock-coupled surface was subtracted.</figcaption></figure>

# 3. Materials and Methods

## 3.1. Protein Production and Purification

GAPA1 coding sequence (excluding the signal peptide) was amplified from a plasmid pda3105 obtained from The RIKEN BioResource Research Center and cloned into pET100/D-TOPO vector (Invitrogen, Waltham, MA, USA)). Insert orientation and sequence was verified by sequencing (Eurofins, Ebersberg, Germany). Recombinant GAPA1 was produced in *E. coli* BL21 (DE3) as a fusion protein with 6xHis-tag and purified by affinity purification with Ni-NTA resin (Thermo Scientific, Waltham, MA, USA). Proteins were buffer-exchanged to phosphate-buffered saline (PBS) by dialysis, concentrated by centrifugation columns (Amicon, Darmstadt, Germany), aliquoted, and stored in 50% (v/v) glycerol at -80 °C for further use.

## 3.2. Protein Mutation

Point mutations to GAPA1 were introduced by replicating the vector using primers that introduce desired substitutions to protein coding sequence (Q5 Site-Directed Mutagenesis Kit, New England Biolabs, Ipswich, MA, USA). Mutations were verified by sequencing (Eurofins, Ebersberg, Germany).

## 3.3. GAPDH Activity Assay

The GAPDH activity of GAPA1 was measured in vitro via colorimetric assay according to the protocol provided by the manufacturer (ab204732, Abcam, Cambridge, UK).

## 3.4. Computational Methods: Starting Structures

For *Arabidopsis thaliana* AtGAPDHA1, we used the crystallographic structure of the photosynthetic homotetrameric (A₄) glyceraldehyde-3-phosphate dehydrogenase A1 (PDB 3K2B) [32]. In silico mutagenesis was performed in CHIMERA [35], using the rotamer tool.
### 3.5. Ab initio Calculations for the Ligand

Molecular geometries of the considered Ligand (SA) have been optimised using the ab initio method at the MP2/6-311++G** level of theory using GAUSSIAN09. Harmonic vibrational frequencies of the different stationary points have been calculated at the same level.

### 3.6. Molecular Dynamics Simulations

Protons were added using CHIMERA [35]. The GROMACS package v5.0.7 [38] was used to run our molecular dynamics simulations. The AMBER99SB-ILDN [39] force field was used to provide molecular mechanics parameters to our proteins. The proteins were put in a cubic cell ("box"), the border of which is at least 1 nm from the protein, and we solvate it with TIP3P [40] explicit water molecules. Counterions ($Na^+$) were added to obtain a neutral system and took the place of water molecules. Then, the protein prepared was first submitted to a quick energy minimisation to ensure that the system has no steric clashes or inappropriate geometry and to avoid unnecessary distortion of the protein when the molecular dynamics simulation is started. The step size of the steepest descent minimisation was 50,000 and minimisation stopped when the maximum force was <1000.0 kJ/mol/nm. Next, we performed two equilibration runs. Bonds will be constrained to enable 2 fs time steps. The first equilibration was performed to reach the target temperature of 300 K in the NVT ensemble (constant Number of particles, Volume, and Temperature) using the modified V-rescale Berendsen thermostat [41] with 100 ps relaxation time. The second equilibration, that is equilibration of pressure, was conducted in the NPT ensemble (wherein the Number of particles, Pressure, and Temperature are all constant) using the Parrinello-Rahman barostat [42] set at 1 bar as reference pressure. The equilibration time was also set at 100 ps relaxation time.

Once our protein was well-equilibrated at the desired temperature and pressure, we release the position restraints and run production MD for data collection. The wild protein (AtGAPDHA1) was subjected to 600 ns simulation with 2 fs time steps. After the simulation was done, we centred the protein in the box and calculate the RMSD (root mean square deviations) and the root mean square fluctuation (RMSF) tools implemented in GROMACS package v5.0.7.

Structures of the protein obtained during simulation were clustered using the Ensemble Cluster tool implemented in the CHIMERA visualisation software [35]. This tool clusters members of a conformational ensemble (i.e., the different structures generated through the simulation) and determines cluster representatives. It is a reimplementation of the method described previously [43].

Docking was performed using AUTODOCK VINA [33,34] implemented as a plug-in in the PYMOL visualisation software. Docking box dimensions were specified so that the protein was surrounded by this box. A flexible docking procedure was applied.

### 3.7. Simulation of the Complex between SA and the Protein

The force field parameters for the ligand were calculated using the ACPYPE software (AnteChamber Python Parser interface) [44]. The charge for SA was *n* = -1. The protein/ligand complex was simulated during at least 20 nsec to check its stability. The complex was built and the unit cell was defined in GROMACS package v5.0.7. After solvation and the addition of ions, the energy of the system was minimised. To do so, first, a position restraint topology was generated for the ligand. Second, the ligand was coupled to the protein for appropriate control of temperature using the make_ndx command of GROMACS package v5.0.7. Equilibration steps and molecular dynamics were done as previously described with energygrps = Protein LIG.

The protein/SA complexes with the highest docking scores were subjected to 20 to 30 ns simulations.

### 3.8. Analysis of the MD Simulations

Visualisation of protein structures with or without ligands were visualised by Chimera [35]. For visualisation reasons, the solvent molecules (water) were not considered.
### 3.9. Surface Plasmon Resonance Analysis

The surface plasmon resonance (SPR) biosensor experiments were performed on a Biacore 3000 instrument (GE Healthcare) of the Platform of Molecular Interactions of the Institute of Biology Paris Seine (IBPS, Sorbonne University). The CM5 sensor chip was used. 3-AESA (called the ligand) was immobilised through primary amino groups to the carboxymethyl dextran matrix of a CM5 sensor chip. 3-AESA (10 mM) diluted in immobilisation buffer (0.1 M Borate) was injected at a flow rate of 10 µL/min for 30 min. The immobilisation was followed by an injection of 70 µL of ethanolamine hydrochloride (1 M) pH 8.5, at a flow rate of 10 µL/min, to saturate the free activated sites of the matrix. A reference surface without ligand was prepared using the same procedure. Buffer HBS-EP (10 mM HEPES, pH 7.4, 150 mM NaCl, 3 mM EDTA, 0.005% surfactant P20) was used as the running buffer and for diluting all the injected molecules. The variations of SPR signal as a function of time (sensorgram) were obtained by passing various concentrations of the injected protein or small molecule (called the analyte) over the ligand surface at a flow rate of 5 µL/min, with an association phase of 5 min and a dissociation phase of 8 min. The surface of the sensor was regenerated with an injection of 25 mM NaOH at a flow rate of 10 µL/min and 30 s contact time. Identical injections over blank surfaces executed in parallel (giving a value of 0 RU) were subtracted from all experiments. The kinetics were evaluated using the BIAevaluation software, version 4.1 (GE Healthcare, Chicago, IL, USA). The data were processed by fitting the binding profiles to a 1:1 Langmuir interaction model with a subsequent extrapolation of response at equilibrium values (Req). The quality of the fit was assessed by the statistical Chi2 value provided by the software (Chi2 values < 10 were considered as acceptable). The fitting of each dataset yielded rates for association ($k_a$) and dissociation ($k_d$), from which the dissociation constant $K_D$ was calculated ($K_D = k_d/k_a$). The $k_{on}$, $k_{off}$, and $K_D$ from 3 experiments were used to calculate the mean values of these variables.

## 4. Discussion and Conclusions

### 4.1. SA Binds Multiple Proteins but No Common Pattern is to Be Expected

SA's ability to bind multiple protein targets and to affect their enzymatic activities is an intriguing phenomenon concerning not only plant physiology but also therapeutic applications of SA derivatives (aspirin). Despite opulent data on the identification of new proteins that bind SA, we still know little about mechanisms of such binding. No conserved SA-binding motif is known, and interaction could occur by various scenarios. For an NRP1, SA binding implicates Cys521/529 and requires the presence of copper ions [9]. No co-crystallisation data of NPR1-SA are yet available. Cys521/529 are not conserved in NPR1 orthologs that nevertheless efficiently bind SA [24]. Not to mention other classes of SABPs. NtSABP2 was co-crystallised with SA and Ser81, Ala13 and His238 residues were producing hydrogen bonds with a carboxylic group of SA [22]. This implies that the characterisation of SA-protein interaction each time will be a bespoke task. Our goal was to propose a new workflow strategy that combines in silico tools with subsequent biochemical validation to study SA-protein interactions at the molecular level.

### 4.2. GAPA1 is a High-Affinity SABP

In this study, we investigated a mode of interaction with SA of a GAPA1 from *Arabidopsis thaliana*, which was previously reported to be an SABP [27]. GAPA1 is a chloroplastic GAPDH that primarily acts in photosynthesis in the Besnon-Calvin cycle.

GAPA1 bound 3-AESA with high affinity $K_D$ = 16.7 nM in SPR analysis. SPR is a powerful approach that provides high-quality data describing the kinetics of protein-ligand interactions as shown for SA and a number of proteins [45,46].

We concluded the interaction to be biologically relevant, since binding competition with SA was much greater than that with 4-HBA, despite the fact that these two molecules are very close from the chemical point of view.
The GAPA1 level of affinity to SA is relatively high when compared to other characterised SA-protein interactions. For example, SABP1 bound SA with a $K_D$ of 14 µM [47], SABP2 bound SA with $K_D$ = 90 nM [21], whereas NPR1 bound SA with $K_D$ = 140 nM [9]. High Mobility Group Box 3 (HMGB3) is an example of an SABP with higher affinity to SA ($K_D$ = 1.5 nM) [46] compared to GAPA1. It has to be mentioned that not the same methods were used to assess the binding of these proteins, so some discrepancy is to be expected.

### 4.3. SA Binds GAPA1 via Arg81 and Asn35

Knowing that SA binds GAPA1, it was important for us to establish which amino acid residues are involved. To do so, we used a protocol that first establishes these residues in molecular dynamics simulations. The first step is protein relaxation in water to obtain representative structures of different clusters; then, each representative ("an average structure") is used for docking calculation. By testing this protocol, we could well replicate the position of SA binding in NtSABP2 where this position is known based on co-crystallisation data [22]. Then, docking poses were assessed by simulating the complex between SA and the protein. In the case of GAPA1, such simulations yielded the identification of a surface pocket that could bind SA. The entrance to this pocket is capped by Arg81 and Thr37. Arg81 was shown to be the residue to which SA establishes a bond in the native protein. This was confirmed by using mutated proteins in the simulations. The importance of Asn35 in the binding was also shown; we suggest that the role of Asn35 is in the structuring of the cavity. By establishing in silico which AA residues could be involved in SA binding, we tested if this is corroborated by in vitro assays. The ability of GAPA1 N35G and R81L to bind 3-AESA in SPR assays (the same protocol as the one to applied WT GAPA1) was indeed significantly reduced.

In this way, our computation approach that use molecular dynamics tools proved to be an efficient strategy for providing primary data as of possible sites of SA-protein interaction. At the same time, we strongly believe that such computational data should always be corroborated by in vitro tests before establishing interaction metrics. In optimistic scenarios, such an approach could greatly speed up finding a place where proteins interact with SA or other ligands for that matter.

### 4.4. SA Inhibits GAPDH Activity of GAPA1. Does the Same Concern Human GAPDH?

In our tests, we have observed that SA inhibits the GAPDH activity of GAPA1 *in vitro*. Inhibition occurred at rather high SA concentrations (2 mM) that are not expected to occur *in planta* (basal SA is around 1 µg/g FW [48]. Still, the biological grounds of this effect are validated by the fact that 4-HBA in the same concentration tested (2 mM) had not induced inhibition. GAPA1 is localised in chloroplasts-organelles that harbour the key pathway for the SA biosynthesis-chorismate pathway. So, one can expect that local SA concentration will be higher in chloroplasts before SA is distributed to cytosol via transporters [49] or passive diffusion. Endogenous SA concentrations are also multifold-increased following SAR [50] or in SA-overaccumulating mutants [51].

Corroborating the fact that SA could indeed be an inhibitor of GAPA1 activity is its localisation in binding complex with the protein. SA is bound exactly at the site of NADP docking the enzyme (Figure 9). In GapB—a GAPDH from spinach—Arg77 (corresponding to Arg81 in GAPA1) was shown to interact with negatively charged 2'-phosphate of NADP [52]. Besides, in GapB R77A, the thioredoxin-mediated redox regulation was abolished and enzyme activity was diminished in vitro [53]. Thus, SA could potentially impede NADP docking to GAPDHs and prevent enzyme catalytic functions/regulation. NADH and SA indeed compete for GAPA1 binding (Figure 11).

The fact that small molecules (such as SA) go into the binding site of a cofactor is not uncommon in enzymology. Cofactor binding sites are targeted by many drugs [54]. For instance, the NAD/NADH cofactor binding site of parasitic S-adenosyl-L-homocysteine hydrolase is targeted for the design of anti-parasitic drugs [55]. The NADH-binding site of the trypanosomal GAPDH was suggested to be a site of choice to develop selective inhibitors against trypanosomes [56]. There are similar concerns
for other cofactors. It is well known that the majority of kinase inhibitors, including dasatinib and PP1, bind at the ATP-binding site [57].

By observing the aligned sequences of GAPDH from different plant species as well as from animals, we have observed that Arg81 and Asn35 are conserved in all GAPDHs (Figure 8). Considering the fact that other plant GAPDH (GAPA2, GAPC1, and GAPC2), as well as human GAPDH also bind SA, it is tempting to speculate that binding occurs at the same position as in GAPA1. This data would be of great importance for analyzing the therapeutic application of SA. In our studies, we have observed that Thr37 could play some role in SA binding (Figure 7). In NAD-dependent GAPDH (either in plants or in animals), position 37 is proline (Figure 8). However, this should not impede potential interactions with SA, since a change of Thr37 to Pro has not led to SA leaving the binding pocket in our simulations (data not shown).

### 4.5. Role of SA Binding and Inhibiting GAPA1 in Plant Stress Responses

However, what could be a biological outcome for SA binding GAPA1 in plants? Chloroplastic GAPDHs are key enzymes of photosynthesis, so SA inhibiting GAPA1 could represent a mechanism for SA to take control of primary metabolism. An observed effect of SA on GAPA1 could similarly have connections to its potential non-canonical activities. For many GAPDH enzymes (mainly glycolytic), numerous so-called moonlighting activities, that are alternative, additional functions in respect to their originally assigned ones, were described [58,59]. GAPC1 is a GAPDH from wheat that directly interacts with the cytochrome b6-f complex iron sulfite subunit (Cytb6f) (in yeast two-hybrid and bimolecular complementation assays). The latter protein has a role in ROS scavenging and could be implicated in drought responses [60].

Moreover, a number of GAPDH enzymes have roles in direct association to pathogen responses. Cytosolic GAPC-1 in *A. thaliana* could directly associate with a negative RNA strand of the tomato bushy stunt virus (TBSV), which was a prerequisite for its asymmetric replication. Binding to SA abolished such activities of GAPC-1 [27]. In *Nicotiana benthamiana*, a chloroplast-localised GAPDH-A interacted with Red clover necrotic mosaic virus (RCNMV) movement protein (MP) [61]. The silencing of NbGAPDH-A resulted in the inhibition of virus multiplication in leaves and prevented MP targeting to the viral replication complex (used by RNA viruses for genome replication) located at a cortical endoplasmic reticulum [61]. A role of SA in this reaction has not been established.

However, no roles outside housekeeping functions were so far reported for GAPA1 from *Arabidopsis thaliana*.

GAPA1-knockout resulted in a slight increase of Arabidopsis resistance to *Pseudomonas syringae* pv. *tomato inoculation* [58], whereas the expression of GAPA1 is significantly reduced following 3 h post-infiltration with 5 µM flg22, which is an elicitor of plant defences [58]. These data provide evidence that GAPA1 could be a negative regulator of immunity. Whether this is due to GAPDH activity or another function of GAPA1, and the role of SA in these reactions, have yet to be established.

**Supplementary Materials:** Supplementary Materials can be found at http://www.mdpi.com/1422-0067/21/13/4678/s1.

**Author Contributions:** All authors have read and agree to the published version of the manuscript. Conceptualisation, I.K. and E.R.; Molecular Dynamics, D.H., I.K. and E.R.; Cloning, D.H. and E.R.; Recombinant protein preparation: D.H., M.H.-M., I.P.; SPR experiments: T.B. and I.P.; Writing—Original Draft Preparation, I.P. and E.R.; Writing—Review & Editing, I.K.; Supervision, E.R.; Funding Acquisition, I.K., E.R. and I.P.; Resources: V.K., L.L.; Methodology: C.E.

**Funding:** Marie Curie/Université Paris-Est Créteil post doc fellowship (Prestige programme) to I.P., PHC DNIPRO grant to E.R. and V.K., M2 grants for D.H. and M.H.-M. from iEES-Paris. OSU Efluve.

**Acknowledgments:** We thank Rahma Dahmani and Halima Mouhib for their help in preparing the ligand PDB file and careful reading of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
References

1. Li, S.; Zhao, J.; Zhai, Y.; Yuan, Q.; Zhang, H.; Wu, X.; Lu, Y.; Peng, J.; Sun, Z.; Lin, L.; et al. The hypersensitive induced reaction 3 (HIR3) gene contributes to plant basal resistance via an EDS1 and salicylic acid-dependent pathway. *Plant J.* **2019**, *98*, 783-797. [CrossRef] [PubMed]

2. Pokotylo, I.; Kravets, V.; Ruelland, E. Salicylic Acid Binding Proteins (SABPs): The Hidden Forefront of Salicylic Acid Signalling. *Int. J. Mol. Sci.* **2019**, *20*, 4377. [CrossRef] [PubMed]

3. Park, S.-W.; Kaimoyo, E.; Kumar, D.; Mosher, S.; Klessig, D.F. Methyl Salicylate Is a Critical Mobile Signal for Plant Systemic Acquired Resistance. *Science* **2007**, *318*, 113-116. [CrossRef] [PubMed]

4. Vlot, A.C.; Liu, P.-P.; Cameron, R.K.; Park, S.-W.; Yang, Y.; Kumar, D.; Zhou, F.; Padukkavidana, T.; Gustafsson, C.; Pichersky, E.; et al. Identification of likely orthologs of tobacco salicylic acid-binding protein 2 and their role in systemic acquired resistance in Arabidopsis thaliana. *Plant J.* **2008**, *56*, 445-456. [CrossRef]

5. Canales, F.J.; Montilla-Bascón, G.; Rispail, N.; Prats, E. Salicylic acid regulates polyamine biosynthesis during drought responses in oat. *Plant Signal. Behav.* **2019**, *14*, e1651183-e1651184. [CrossRef]

6. Jayakannan, M.; Bose, J.; Babourina, O.; Rengel, Z.; Shabala, S. Salicylic acid improves salinity tolerance in Arabidopsis by restoring membrane potential and preventing salt-induced K+ loss via a GORK channel. *J. Exp. Bot.* **2013**, *64*, 2255-2268. [CrossRef]

7. Needs, C.J.; Brooks, P.M.; Brooks, P.P.M. Clinical Pharmacokinetics of the Salicylates. *Clin. Pharmacokinet.* **1985**, *10*, 164-177. [CrossRef]

8. Janda, M.; Ruelland, E. Magical mystery tour: Salicylic acid signalling. *Environ. Exp. Bot.* **2015**, *114*, 117-128. [CrossRef]

9. Wu, Y.; Zhang, D.; Chu, J.Y.; Boyle, P.; Wang, Y.; Brindle, I.; De Luca, V.; Després, C. The Arabidopsis NPR1 Protein Is a Receptor for the Plant Defense Hormone Salicylic Acid. *Cell Rep.* **2012**, *1*, 639-647. [CrossRef]

10. Glazebrook, J.; Rogers, E.E.; Ausubel, F.M. Isolation of Arabidopsis Mutants with Enhanced Disease Susceptibility by Direct Screening. *Genetics* **1996**, *143*, 973-982.

11. Blanco, F.; Garreton, V.; Frey, N.; Domínguez, C.; Perez-Acle, T.; Van Der Straeten, D.; Jordana, X.; Holuigue, L.; Blanco-Herrera, F. Identification of NPR1-Dependent and Independent Genes Early Induced by Salicylic Acid Treatment in Arabidopsis. *Plant Mol. Boil.* **2005**, *59*, 927-944. [CrossRef] [PubMed]

12. Canet, J.V.; Dobon, A.; Roig, A.; Tornero, P. Structure-function analysis of npr1 alleles in Arabidopsis reveals a role for its paralogs in the perception of salicylic acid. *Plant Cell Environ.* **2010**, *33*, 1911-1922. [CrossRef]

13. Rong, D.; Luo, N.; Mollet, J.-C.; Liu, X.; Yang, Z. Salicylic Acid Regulates Pollen Tip Growth through an NPR3/NPR4-Independent Pathway. *Mol. Plant* **2016**, *9*, 1478-1491. [CrossRef]

14. Tan, S.; Abas, M.; Verstraeten, I.; Glanc, M.; Molnár, G.; Hajný, J.; Lasák, P.; Petřík, I.; Russinova, E.; Petrášek, J.; et al. Salicylic Acid Targets Protein Phosphatase 2A to Attenuate Growth in Plants. *Curr. Boil.* **2020**, *30*, 381-395. [CrossRef] [PubMed]

15. Lai, Y.-S.; Renna, L.; Yarema, J.; Ruberti, C.; He, S.Y.; Brandizzi, F. Salicylic acid-independent role of NPR1 is required for protection from proteotoxic stress in the plant endoplasmic reticulum. *Proc. Natl. Acad. Sci. USA* **2018**, *115*, E5203-E5212. [CrossRef] [PubMed]

16. Mou, Z.; Fan, W.; Dong, X. Inducers of Plant Systemic Acquired Resistance Regulate NPR1 Function through Redox Changes. *Cell* **2003**, *113*, 935-944. [CrossRef]

17. Tada, Y.; Spoel, S.H.; Pajerowska-Mukhtar, K.; Mou, Z.; Song, J.; Wang, C.; Zuo, J.; Dong, X. Plant Immunity Requires Conformational Charges of NPR1 via S-Nitrosylation and Thioredoxins. *Science* **2008**, *321*, 952-956. [CrossRef] [PubMed]

18. Krinke, O.; Ruelland, E.; Valentová, O.; Vergnolle, C.; Renou, J.-P.; Taconnat, L.; Flemr, M.; Burketová, L.; Zachowski, A. Phosphatidylinositol 4-Kinase Activation Is an Early Response to Salicylic Acid in Arabidopsis Suspension Cells. *Plant Physiol.* **2007**, *144*, 1347-1359. [CrossRef] [PubMed]

19. Krinke, O.; Flemr, M.; Vergnolle, C.; Collin, S.; Renou, J.-P.; Taconnat, L.; Yu, A.; Burketová, L.; Valentová, O.; Zachowski, A.; et al. Phospholipase D Activation Is an Early Component of the Salicylic Acid Signaling Pathway in Arabidopsis Cell Suspensions. *Plant Physiol.* **2009**, *150*, 424-436. [CrossRef]

20. Chen, Z.; Ricigliano, J.W.; Klessig, D.F. Purification and characterization of a soluble salicylic acid-binding protein from tobacco. *Proc. Natl. Acad. Sci. USA* **1993**, *90*, 9533-9537. [CrossRef]

21. Du, H.; Klessig, D.F. Identification of a Soluble, High-Affinity Salicylic Acid-Binding Protein in Tobacco. *Plant Physiol.* **1997**, *113*, 1319-1327. [CrossRef] [PubMed]
22. Forouhar, F.; Yang, Y.; Kumar, D.; Chen, Y.; Fridman, E.; Park, S.W.; Chiang, Y.; Acton, T.B.; Montelione, G.T.; Pichersky, E.; et al. Structural and biochemical studies identify tobacco SABP2 as a methyl salicylate esterase and implicate it in plant innate immunity. *Proc. Natl. Acad. Sci. USA* **2005**, *102*, 1773–1778. [CrossRef] [PubMed]

23. Slaymaker, D.H.; Navarre, D.A.; Clark, D.; Del Pozo, O.; Martin, G.B.; Klessig, D.F. The tobacco salicylic acid-binding protein 3 (SABP3) is the chloroplast carbonic anhydrase, which exhibits antioxidant activity and plays a role in the hypersensitive defense response. *Proc. Natl. Acad. Sci. USA* **2002**, *99*, 11640–11645. [CrossRef] [PubMed]

24. Fu, Z.Q.; Yan, S.; Saleh, A.; Wang, W.; Ruble, J.; Oka, N.; Mohan, R.; Spoel, S.H.; Tada, Y.; Zheng, N.; et al. NPR3 and NPR4 are receptors for the immune signal salicylic acid in plants. *Nature* **2012**, *486*, 228–232. [CrossRef]

25. Medina-Puche, L.; Castelló, M.J.; Canet, J.V.; Lamilla, J.; Colombo, M.L.; Tornero, P. β-carbonic anhydrases play a role in salicylic acid perception in Arabidopsis. *PLoS ONE* **2017**, *12*, e0181820. [CrossRef] [PubMed]

26. Manohar, M.; Tian, M.; Moreau, M.; Park, S.-W.; Choi, H.W.; Fei, Z.; Friso, G.; Asif, M.; Manosalva, P.; Von Dahl, C.C.; et al. Identification of multiple salicylic acid-binding proteins using two high throughput screens. *Front. Plant Sci.* **2015**, *5*, 5. [CrossRef]

27. Tian, M.; Sasvári, Z.; Gonzalez, P.A.; Friso, G.; Rowland, E.; Liu, X.-M.; Van Wijk, K.J.; Nagy, P.D.; Klessig, D.F.; Sasvani, Z. Salicylic Acid Inhibits the Replication of Tomato bushy stunt virus by Directly Targeting a Host Component in the Replication Complex. *Mol. Plant-Microbe Interact.* **2015**, *28*, 379–386. [CrossRef] [PubMed]

28. Dachineni, R.; Ai, G.; Kumar, D.R.; Sadhu, S.S.; Tummala, H.; Bhat, G. Cyclin A2 and CDK2 as Novel Targets of Aspirin and Salicylic Acid: A Potential Role in Cancer Prevention. *Mol. Cancer Res.* **2015**, *14*, 241–252. [CrossRef]

29. Choi, H.W.; Tian, M.; Manohar, M.; Harraz, M.M.; Park, S.-W.; Schroeder, F.C.; Snyder, S.H.; Klessig, D.F. Human GAPDH Is a Target of Aspirin's Primary Metabolite Salicylic Acid and Its Derivatives. *PLoS ONE* **2015**, *10*, e0143447. [CrossRef]

30. Chen, Z.; Silva, H.; Klessig, D. Active oxygen species in the induction of plant systemic acquired resistance by salicylic acid. *Science* **1993**, *262*, 1883–1886. [CrossRef]

31. Moreau, M.; Westlake, T.; Zampogna, G.; Popescu, G.V.; Tian, M.; Noutsos, C.; Popescu, S. The Arabidopsis oligopeptidases TOP1 and TOP2 are salicylic acid targets that modulate SA-mediated signaling and the immune response. *Plant J.* **2013**, *76*, 603–614. [CrossRef] [PubMed]

32. Fermani, S.; Sparla, F.; Marri, L.; Thumiger, A.; Pupillo, P.; Falini, G.; Trost, P. Structure of photosynthetic glyceraldehyde-3-phosphate dehydrogenase (isoform A4) from Arabidopsis thaliana in complex with NAD. *Acta Crystallogr. Sect. F Struct. Boil. Cryst. Commun.* **2010**, *66*, 621–626. [CrossRef] [PubMed]

33. Seeliger, D.; De Groot, B.L. Ligand docking and binding site analysis with PyMOL and Autodock/Vina. *J. Comput. Mol. Des.* **2010**, *24*, 417–422. [CrossRef] [PubMed]

34. Trott, O.; Olson, A.J. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. *J. Comput. Chem.* **2010**, *31*, 455–461. [CrossRef]

35. Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Couch, G.S.; Greenblatt, D.M.; Meng, E.C.; Ferrin, T.E. UCSF Chimera—A visualization system for exploratory research and analysis. *J. Comput. Chem.* **2004**, *25*, 1605–1612. [CrossRef]

36. Laskowski, R.A.; Swindells, M.B. LigPlot+: Multiple Ligand–Protein Interaction Diagrams for Drug Discovery. *J. Chem. Inf. Model.* **2011**, *51*, 2778–2786. [CrossRef]

37. Robert, X.; Gouet, P. Deciphering key features in protein structures with the new ENDscript server. *Nucleic Acids Res.* **2014**, *42*, W320–W324. [CrossRef] [PubMed]

38. Abraham, M.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J.C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. *SoftwareX* **2015**, *1*, 19–25. [CrossRef]

39. Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J.L.; Dror, R.O.; Shaw, D.E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. *Proteins: Struct. Funct. Bioinform.* **2010**, *78*, 1950–1958. [CrossRef] [PubMed]

40. Jorgensen, W.L. Quantum and statistical mechanical studies of liquids. 10. Transferable intermolecular potential functions for water, alcohols, and ethers. Application to liquid water. *J. Am. Chem. Soc.* **1981**, *103*, 335–340. [CrossRef]
41. Berendsen, H.J.C.; Postma, J.P.M.; DiNola, A.; Haak, J.R.; Van Gunsteren, W.F. Molecular dynamics with coupling to an external bath. *J. Chem. Phys.* **1984**, *81*, 3684. [CrossRef]

42. Parrinello, M. Polymorphic transitions in single crystals: A new molecular dynamics method. *J. Appl. Phys.* **1981**, *52*, 7182. [CrossRef]

43. Kelley, L.A.; Gardner, S.P.; Sutcliffe, M. An automated approach for clustering an ensemble of NMR-derived protein structures into conformationally related subfamilies. *Protein Eng. Des. Sel.* **1996**, *9*, 1063–1065. [CrossRef] [PubMed]

44. Da Silva, A.W.S.; Vranken, W. ACPYPE - AnteChamber Python Parser interfacE. *BMC Res. Notes* **2012**, *5*, 367. [CrossRef] [PubMed]

45. Manohar, M.; Choi, H.W.; Manosalva, P.; Austin, C.A.; Peters, J.E.; Klessig, D.F. Plant and Human MORC Proteins Have DNA-Modifying Activities Similar to Type II Topoisomerases, but Require One or More Additional Factors for Full Activity. *Mol. Plant-Microbe Interact.* **2017**, *30*, 87–100. [CrossRef] [PubMed]

46. Choi, H.W.; Manohar, M.; Manosalva, P.; Tian, M.; Moreau, M.; Klessig, D.F. Activation of Plant Innate Immunity by Extracellular High Mobility Group Box 3 and Its Inhibition by Salicylic Acid. *PLoS Pathog.* **2016**, *12*, e1005518. [CrossRef] [PubMed]

47. Chen, Z.; Klessig, D.F. Identification of a soluble salicylic acid-binding protein that may function in signal transduction in the plant disease-resistance response. *Proc. Natl. Acad. Sci. USA* **1991**, *88*, 8179–8183. [CrossRef] [PubMed]

48. Dempsey, D.A.; Klessig, D.F. How does the multifaceted plant hormone salicylic acid combat disease in plants and are similar mechanisms utilized in humans? *BMC Boil.* **2017**, *15*, 23. [CrossRef]

49. Serrano, M.; Wang, B.; Aryal, B.; Garcion, C.; Mansour, E.A.; Heck, S.; Geisler, M.; Mauch, F.; Nawrath, C.; Metraux, J.P. Export of salicylic acid from the chloroplast requires the multidrug and toxin extrusion-like transporter EDS5. *Plant Physiol.* **2013**, *162*, 1815–1821. [CrossRef]

50. Wildermuth, M.C.; Dewdney, J.; Wu, G.; Ausubel, F.M. Isochorismate synthase is required to synthesize salicylic acid for plant defence. *Nature* **2001**, *414*, 562–565. [CrossRef]

51. Šašek, V.; Janda, M.; Delage, E.; Puyaubert, J.; Guivarc'H, A.; Maseda, E.L.; Dobrev, P.; Caius, J.; Bóka, K.; Valentová, O.; et al. Constitutive salicylic acid accumulation in ppi4kIIIβ1β2Arabidopsis plants stunts rosette but not root growth. *New Phytol.* **2014**, *203*, 805–816. [CrossRef] [PubMed]

52. Sparla, F.; Fermani, S.; Falini, G.; Zaffagnini, M.; Ripamonti, A.; Sabatino, P.; Pupillo, P.; Trost, P. Coenzyme Site-directed Mutants of Photosynthetic A4-GAPDH Show Selectively Reduced NADPH-dependent Catalysis, Similar to Regulatory AB-GAPDH Inhibited by Oxidized Thioredoxin. *J. Mol. Boil.* **2004**, *340*, 1025–1037. [CrossRef]

53. Sparla, F.; Zaffagnini, M.; Wedel, N.; Scheibe, R.; Pupillo, P.; Trost, P. Regulation of Photosynthetic GAPDH Dissected by Mutants1. *Plant Physiol.* **2005**, *138*, 2210–2219. [CrossRef]

54. Hua, Y.H.; Wu, C.Y.; Sargsyan, K.; Lim, C. Sequence-motif Detection of NAD(P)-binding Proteins: Discovery of a Unique Antibacterial Drug Target. *Sci. Rep.* **2014**, *4*, 6471. [CrossRef]

55. Cai, S.; Li, Q.-S.; Fang, J.; Borchardt, R.T.; Kuczera, K.; Middaugh, C.R.; Schowen, R.L. The rationale for targeting the NAD/NADH cofactor binding site of parasitic S-adenosyl-L-homocysteine hydrolase for the design of anti-parasitic drugs. *Nucleosides Nucleotides Nucleic Acids* **2009**, *28*, 485–503. [CrossRef] [PubMed]

56. Suresh, S.; Turley, S.; Opperdoes, F.; Michels, P.A.; Hol, W.G.J. A potential target enzyme for trypanocidal drugs revealed by the crystal structure of NAD-dependent glycerol-3-phosphate dehydrogenase from Leishmania mexicana. *Structure* **2000**, *8*, 541–552. [CrossRef]

57. Shan, Y.; Kim, E.T.; Eastwood, M.P.; Dror, R.O.; Seeliger, M.; Shaw, D.E. How Does a Drug Molecule Find Its Target Binding Site? *J. Am. Chem. Soc.* **2011**, *133*, 9181–9183. [CrossRef] [PubMed]

58. Henry, E.; Fung, N.; Liu, J.; Drakakaki, G.; Coaker, G. Beyond Glycolysis: GAPDHs Are Multi-functional Enzymes Involved in Regulation of ROS, Autophagy, and Plant Immune Responses. *PLoS Genet.* **2015**, *11*, e1005199. [CrossRef]

59. Zaffagnini, M.; Fermani, S.; Costa, A.; Lemaire, S.D.; Trost, P. Plant cytoplasmic GAPDH: Redox post-translational modifications and moonlighting properties. *Front. Plant Sci.* **2013**, *4*, 450. [CrossRef]
60. Li, X.; Wei, W.; Li, F.; Zhang, L.; Deng, X.; Liu, Y.; Yang, S. The Plastidial Glyceraldehyde-3-Phosphate Dehydrogenase Is Critical for Abiotic Stress Response in Wheat. *Int. J. Mol. Sci.* **2019**, *20*, 1104. [CrossRef]

61. Kaido, M.; Abe, K.; Mine, A.; Hyodo, K.; Taniguchi, T.; Taniguchi, H.; Mise, K.; Okuno, T. GAPDH-A Recruits a Plant Virus Movement Protein to Cortical Virus Replication Complexes to Facilitate Viral Cell-to-Cell Movement. *PLoS Pathog.* **2014**, *10*, e1004505. [CrossRef] [PubMed]

<figure><img src="image_15.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).