Respiratory and C₄-photosynthetic NAD-malic enzyme coexist in bundle sheath cell mitochondria and evolved via association of differentially adapted subunits

Abstract

In plant mitochondria, nicotinamide adenine dinucleotide-malic enzyme (NAD-ME) has a housekeeping function in malate respiration. In different plant lineages, NAD-ME was independently co-opted in C₄ photosynthesis. In the C₄ Cleome species, Gynandropsis gynandra and Cleome angustifolia, all NAD-ME genes (NAD-MEα, NAD-MEβ1, and NAD-MEβ2) were affected by C₄ evolution and are expressed at higher levels than their orthologs in the C₃ species Tarenaya hassleriana. In T. hassleriana, the NAD-ME housekeeping function is performed by two heteromers, NAD-MEα/β1 and NAD-MEα/β2, with similar biochemical properties. In both C₄ species, this role is restricted to NAD-MEα/β2. In the C₄ species, NAD-MEα/β1 is exclusively present in the leaves, where it accounts for most of the enzymatic activity. Gynandropsis gynandra NAD-MEα/β1 (GgNAD-MEα/β1) exhibits high catalytic efficiency and is differentially activated by the C₄ intermediate aspartate, confirming its role as the C₄-decarboxylase. During C₄ evolution, NAD-MEβ1 lost its catalytic activity; its contribution to the enzymatic activity results from a stabilizing effect on the associated α-subunit and the acquisition of regulatory properties. We conclude that in bundle sheath cell mitochondria of C₄ species, the functions of NAD-ME as C₄ photosynthetic decarboxylase and as a housekeeping enzyme coexist and are performed by isoforms that combine the same α-subunit with differentially adapted β-subunits.

Introduction

Many of the world’s most productive crops are species that exhibit C₄ photosynthetic metabolism, as an adaptation to high light intensities, temperatures, and dryness (Sage et al., 2012). C₄ plants have evolved biochemical pumps that concentrate CO₂ at the site of ribulose-1,5-bisphosphate carboxylase/oxygenase, which decrease the oxygenation reaction and, thus, limits the wasteful flux through photorespiration (Furbank and Hatch, 1987). Compared with plants using the ancestral C₃ photosynthetic pathway, C₄ plants show higher water and nitrogen use efficiency, which allows increased productivity in warm habitats (Edwards et al., 2010). The transition from C₃ to C₄ photosynthesis involved complex alterations to leaf anatomy and biochemistry. This occurred in at least 66 linages of angiosperms independently and, thus, is an example of convergent evolution (Sage, 2016). Most C₄ plants developed a spatial separation of the biochemical components of the CO₂ pump by utilizing two adjacent cell types. In these plants, phosphoenolpyruvate (PEP) carboxylase is located in the cytosol of mesophyll cells, and the formed C₄ intermediates are shuttled to the bundle sheath cells (BSCs) (Drincovich et al., 2011). The release of CO₂ from malate in BSCs is mediated mainly by two different decarboxylases: nicotinamide adenine dinucleotide phosphate (NADP)-malic enzyme (NADP-ME, EC 1.1.1.40), which is exclusively located to chloroplasts, and NAD-ME (EC 1.1.1.39), with exclusive location in the mitochondria (Maier et al., 2011; Wang et al., 2014). The majority of the C₄ linages mainly use NADP-ME (Sage et al., 2011). The use of NAD-ME is most widespread in eudicot species, where it is found in approximately 20 C₄ lineages (Sage et al., 2011).

The C₄ enzymes did not evolve de novo. Instead, they were recruited from existing housekeeping isoforms (Monson, 2003; Aubry et al., 2011; Maier et al., 2011; Christin et al., 2013a, 2013b). Most C₄ enzymes originated by gene duplication of an existing nonphotosynthetic isoform and the subsequent neofunctionalization of one gene copy (Tausta et al., 2002; Saigo et al., 2004; Christin et al., 2013a, 2013b; Ludwig, 2016a, 2016b; Alvarez et al., 2019; Bovdilova et al., 2019). The C₄-NAD-ME turned out to be an exceptional case. Recently, we described that the evolution of NAD-ME in land plants is marked by sub-functionalization and differences in the frequency of gene duplication of the two paralogous α- and β-NAD-ME gene lineages (Tronconi et al., 2020). Most angiosperm genomes maintained a 1:1 α-NAD-ME/β-NAD-ME relative gene dosage, but a significantly high proportion of species with C₄-NAD-ME-type photosynthesis have a non-1:1 ratio of α-NAD-ME/β-NAD-ME. Specifically, some C₃ and C₄ species of Brassicales possess a single NAD-MEα gene and two NAD-MEβ genes (β1 and β2) (Tronconi et al., 2020). This is the case of the two independently evolved C₄ species of the genus Cleome, Gynandropsis gynandra and Cleome angustifolia, which not only developed different types of the Kranz anatomy (Koteyeva et al., 2011) but also originated in different geographic areas (Feodorova et al., 2010). In these C₄ species, all three NAD-ME genes were affected by C₄ evolution with the encoded NAD-MEβ1 subunits exhibiting several amino acids identically substituted and positively selected in both C₄ species (Tronconi et al., 2020). Likely, these changes provided the basis for the recruitment of NAD-ME and its function in C₄-physiology.

In its housekeeping role, plant NAD-ME is likely involved in mitochondrial malate respiration, functioning as an associated enzyme to the tricarboxylic acid (TCA) cycle (Grover et al., 1981; Artus and Edwards, 1985; Tronconi et al., 2008; Fuchs et al., 2020). Plant NAD-ME mostly operates as a heteromer of α-(∼63 kDa) and β- (∼58 kDa) subunits that share about 65% sequence identity (Tronconi et al., 2008), which are encoded by the paralogous α- and β-NAD-ME genes (Tronconi et al., 2020). Depending on the species and/or plant tissue, NAD-ME is active as a dimer, tetramer, and octamer (Artus and Edwards, 1985). In the C₃ species Arabidopsis thaliana, NAD-ME is a heterodimer of the α-(NAD-ME1, At4g13560) and β-(NAD-ME2, At4g00570) subunits, while a small proportion of homodimers is present in components of the floral organ (Maier et al., 2011; Tronconi et al., 2018). The recombinant homodimers and the reconstituted heterodimer present similar catalytic efficiencies but differ in kinetic mechanisms and their regulation by metabolic effectors (Tronconi et al., 2008, 2010a, 2010b, 2015). It is intriguing that there are no references to housekeeping or photosynthetic isoforms of NAD-ME in any C₄ plants. Due to the intricate molecular mechanism of evolution of NAD-ME (Tronconi et al., 2020), it is not surprising that despite its central importance in the C₄-pathway, C₄-NAD-ME has never been characterized at the molecular level. Even comparative transcriptomes of leaves of the closely related Cleome Tarenaya hassleriana, a C₃ ornamental species, and G. gynandra, an NAD-ME C₄ orphan crop species (Sogbohossou et al., 2018), did not reveal any C₄-specific NAD-ME transcript (Brautigam et al., 2011).

Elucidation of the structural composition and biochemical characteristics of the components of the C₄ pathway is not only important to gain knowledge on the evolutionary mechanisms underlying their origin but also because many efforts are currently underway to install characteristics of C₄ photosynthesis in leaves of C₃ crops (Kajala et al., 2011; Ermakova et al., 2020; Lin et al., 2020). This endeavor can only be fruitful if it is built on a detailed understanding of the particular characteristics of the components of the C₄ pathway. In this work, we addressed how NAD-ME has been adapted to perform the housekeeping and C₄ functions in C₄ plant mitochondria in two independently evolved C₄ species of the genus Cleome, G. gynandra and C. angustifolia. We found that in C₄ Cleome, the functions of NAD-ME as a respiratory enzyme and specific photosynthetic decarboxylase are performed by separated enzymatic entities originating from the combination of the α-subunit with differentially evolved β-subunits.

Results

Expression of NAD-ME genes in C₃ and C₄ Cleome species

To find candidate gene(s) of C₄-NAD-ME in Cleome, we analyzed the diel expression pattern of all NAD-ME genes (α, β1, and β2) in fully expanded leaves of T. hassleriana (C₃) and G. gynandra (C₄) by quantitative transcriptional analysis (Figure 1;Supplemental Table S1) using the respective ACT2 gene as reference.

Figure 1

Expression levels of genes encoding NAD-ME proteins in G. gynandra (pink) and T. hassleriana (blue) relative to ACT2 in leaves of 5-week-old plants over a full day. The values represent the mean ± se (standard error); n=three independent replicates. Within species, the same letters indicate no statistical differences, different letters indicate statistical differences among the time points based on two-way analysis of variance (ANOVA) for gene and time as factors (α =  0.05) with a post-hoc Tukey’s test (results of multiple comparisons with confidence intervals and significance levels are shown in Supplemental Table S1).

In both species, the expression of the NAD-MEα genes slightly peaks in the middle of the light period (Figure 1). The expression patterns of the NAD-MEβ genes do not change drastically during a diurnal cycle in the C₃ and C₄ species (Figure 1). We found that GgNAD-MEα is expressed 2- to 10-fold higher than GgACT2, while GgNAD-MEβ1 reaches up to 3.7-fold and GgNAD-MEβ2 up to 1.3-fold higher expression than GgACT2 at midday (Figure 1). ThNAD-MEα expression levels are similar to that of ThACT2, while ThNAD-MEβ1 and ThNAD-MEβ2 show approximately 20-fold lower expression than the control gene (Figure 1). As all three GgNAD-ME genes are expressed at much higher levels than their homologs in T. hassleriana, we cannot unambiguously identify a concrete C₄-NAD-ME candidate gene by the transcriptional analysis.

NAD-ME in photosynthetic and heterotrophic organs of C₃ and C₄ Cleome

We explored the presence of NAD-ME isoforms in photosynthetic and heterotrophic organs of C₃ (T. hassleriana) and C₄ (G. gynandra and C. angustifolia) Cleome species by coupling in gel NAD-ME activity assay after gradient native polyacrylamide gel electrophoresis (PAGE) with the identification of NAD-ME subunits in the gel slices through mass spectroscopy (MS). Overall, in gel activity assays indicated different protein bands showing NAD-ME activity distributed over species and different organs (migration levels 1–4 in Figure 2).

Figure 2

Behavior of NAD-ME entities of Cleome species. A, Gradient native PAGE coupled to NAD-ME activity assay of soluble protein extracts of leaves (L, 40 µg protein) of C. angustifolia (Ca) and G. gynandra (Gg). Representative gel of n = 4. B, Gradient native PAGE coupled to NAD-ME activity assay of soluble protein extracts of leaves (L, 40 µg protein) and roots (R, 20 µg protein) of G. gynandra (Gg), and T. hassleriana (Th); A. thaliana (At) soluble protein extracts of leaves were used as control (40-µg protein) (Tronconi et al., 2008). For the in gel activity assays in A and B, the gels were incubated for 150 min in the assay medium at pH 6.8. A violet precipitate indicates NAD-ME activity. Migration levels 1–4 are indicated between A and B. The molecular weight of the control AtNAD-ME (Tronconi et al., 2008) is indicated on the right. Representative gel of n = 10. C, Kinetic parameters of recombinant NAD-ME entities identified in protein extracts of G. gynandra and T. hassleriana. Kinetic data were fitted by nonlinear regression. ^aData taken from (Tronconi et al., 2010). Values represent mean ± se; n=three independent enzyme preparations, each measured in triplicate. As GgNAD-MEβ1 is not active, k_cat was calculated assuming the formation of a heterodimer with only one active site in GgNAD-MEα/β1.

Leaves of G. gynandra and C. angustifolia show two main protein bands with NAD-ME activity of very different mobilities (migration levels 1 and 4, Figure 2A); these bands are not present in the C₃ Cleome species (Figure 2B). One protein band with NAD-ME activity has a very high mobility (migration level 4, Figure 2A) and showed the presence of only NAD-MEα and NAD-MEβ1 after MS analysis (Supplemental Table S2). As the mobility of this protein band is similar to that of the heterodimer NAD-MEα/β found in A. thaliana leaves (Tronconi et al., 2008; Figure 2B), we propose that it corresponds to a NAD-MEα/β1 heterodimer. The second protein band found exclusively in leaves of C₄ species has very low mobility (migration level 1, Figure 2A) and appears rapidly after incubation in the activity assay buffer (Supplemental Figure S1A). MS analysis of this low mobility protein band from both C₄ species showed the presence of all three NAD-ME subunits (Supplemental Tables S2). Most likely, two different heteromers with a similar velocity of migration, NAD-MEα/ß1, and NAD-MEα/ß2, contribute to the NAD-ME activity of this low mobility band of leaves of the C₄ species.

We further analyzed roots of the C₄ species G. gynandra, which show a predominant NAD-ME activity band with low mobility at migration level 1 (Figure 2B). MS analysis of this protein band identified peptides corresponding to NAD-MEα and NAD-MEβ2 (Supplemental Table S2), indicating that the NAD-MEα/β2 heteromer is the only isoform present in roots.

Leaves and roots of the C₃ plant T. hassleriana show a protein band with NAD-ME activity that is not found in organs of the C₄ species (migration level 3, Figure 2B). This protein band has a slightly lower mobility than that of the A. thaliana NAD-ME heterodimer. MS analysis of the protein bands indicated the presence of NAD-MEα and identified peptides that correspond to both NAD-ME β-subunits (Supplemental Table S2). Thus, we suggest that the NAD-MEα/β1 and NAD-MEα/β2 isoforms are present in the protein bands of leaves and roots and have very similar electrophoretic mobilities.

The three analyzed species, G. gynandra, C. angustifolia, and T. hassleriana, share an additional protein band with NAD-ME activity and intermediate mobility (migration level 2, Figure 2, A and B). MS analysis indicated the presence of a NADP-ME in all these protein bands, which was confirmed by an in gel NADP-ME activity assay after native PAGE (Supplemental Figure S1B). As the identified NADP-MEs (Th09755, Th20644, and Gg17666) are homologs of co-factor promiscuous cytosolic AtNADP-MEs (Gerrard Wheeler et al., 2008), we conclude that the band at migration level 2 found in most samples represents a NADP-ME that also uses NAD as a cofactor at least in vitro.

Taken together, leaves and roots of the Cleome C₃ species T. hassleriana most likely possess two NAD-ME isoforms, ThNAD-MEα/β1 and ThNAD-MEα/β2, formed by the alternative combination of the α-subunit with either of the β-subunits. The Cleome C₄ species G. gynandra and C. angustifolia also possess two NAD-ME isoforms, Ca/GgNAD-MEα/β1, and Ca/GgNAD-MEα/β2. Ca/GgNAD-MEα/β2 is present in heterotrophic and photosynthetic organs and is assembled as a heteromer of low mobility. Ca/GgNAD-MEα/β1 is exclusively present in photosynthetic tissues. In the conditions of our assay, Ca/GgNAD-MEα/β1 is responsible for the major part of NAD-ME activity and is active in two structural assemblies: as a heterodimer and as a higher-order heteromer.

Biochemical properties of recombinant NAD-ME of C₃ and C₄ Cleome

To analyze if the identified NAD-ME entities in C₄ Cleome differ in biochemical properties, we performed in vitro analyses with GgNAD-MEs. We co-expressed the mature GgNAD-ME subunits in Escherichia coli and purified them using the His-Tag fused to only one subunit (Supplemental Figure S2). That way, the co-elution of the subunits indicates specific interactions of the proteins. We found that GgNAD-MEα interacts with the NAD-MEβ subunits to form the GgNAD-MEα/β1 and GgNAD-MEα/β2 heteromers, while the recombinant GgNAD-MEβ1 and GgNAD-MEβ2 subunits do not interact with each other (Supplemental Figure S2). These results reinforced the results of the in gel NAD-ME assay coupled to MS analysis, showing that in leaves of C₄ Cleome, NAD-ME forms two heteromers constituted by the association of an α-subunit with a different β-subunit (GgNAD-MEα/β1 and GgNAD-MEα/β2, Figure 2A).

The recombinant GgNAD-MEα/β1 and GgNAD-MEα/β2 are active enzymes with a similar optimum pH of 6.2–6.4 and K_m,NAD of 0.6–0.7 mM, but different kinetic parameters with respect to the substrate malate. GgNAD-MEα/β2 has kinetic parameters in the same range as that of the Arabidopsis isoform (Figure 2C). In contrast, GgNAD-MEα/β1 has an approximately 30 times higher catalytic efficiency than GgNAD-MEα/β2 (Figure 2C). This high catalytic efficiency results from a high apparent affinity to malate, which is 12 times greater than that of GgNAD-MEα/β2, and a turnover rate 2 times higher than that of GgNAD-MEα/β2 (Figure 2C). Similarly, we also expressed and analyzed the NAD-ME isoforms identified in the C₃ species T. hassleriana (ThNAD-MEα/β1 and ThNAD-MEα/β2). We found that both isoforms have similar kinetic parameters to those of GgNAD-MEα/β2 and AtNAD-MEα/β (Figure 2C), indicating their housekeeping functions in mitochondrial respiration.

In gradient native PAGE, the recombinant GgNAD-MEα/β1 is found as two bands of different electrophoretic mobilities, which let us assume that in the conditions of the assay, this isoform can be found as a heterodimer and heterotetramer (Supplemental Figure S3), as also observed with extracts of Cleome organs (Figure 2A). Consistently, GgNAD-MEα/β1 eluted into two peaks after gel filtration chromatography: a peak representing the dimer (molecular mass of 132 ± 15 kDa) and a second broad peak representing proteins of high molecular weight (between 250 and 460 kDa). GgNAD-MEα/β1 recovered in both protein peaks presented NAD-ME activity, with the specific activity of the high molecular weight protein assembly being almost two times higher than that of the heterodimer. Recombinant GgNAD-MEα/β2 eluted as a broad peak, consistent with a high molecular weight protein assembly (between 250 and 460 kDa).

Taken together, the analysis of recombinant NAD-ME proteins indicates that in G. gynandra, the association of the α- and β1-subunits renders a NAD-ME isoform, GgNAD-MEα/β1, with a high affinity for malate and high catalytic efficiency. These are expected properties for an enzyme to fulfill a role in a high-flux metabolic pathway, such as the C₄ photosynthetic pathway. On the other hand, the housekeeping NAD-ME of G. gynandra is formed by the association of α- and β2-subunits, GgNAD-MEα/β2. The C₃ species T. hassleriana possesses two housekeeping NAD-ME isoforms with similar biochemical properties, which are formed by the alternative combination of the α-subunit with either of the β-subunits, ThNAD-MEα/β1 and ThNAD-MEα/β2.

NAD-MEβ1 evolved as a non-catalytic subunit in C₄ Cleome

To aid in understanding the molecular determinants of the evolution of different NAD-ME entities in Cleome, we analyzed the single NAD-ME proteins of G. gynandra and T. hassleriana produced in vitro (Figure 3A;Supplemental Figure S4). The results indicated that all recombinant proteins, except for one, are active enzymes (Figure 3A). Intriguingly, we found that GgNAD-MEβ1 is the only inactive protein. Similarly, we found that the recombinant NAD-MEβ1 of C. angustifolia is also inactive.

Figure 3

NAD-MEβ1 evolved as a noncatalytic subunit in C₄ Cleome. A, Kinetic parameters of recombinant NAD-ME single proteins from G. gynandra, C. angustifolia, and T. hassleriana. Kinetic data were fitted by nonlinear regression. Values represent mean ± se of at least three independent enzyme preparations. ^aData taken from (Tronconi et al., 2008). B, Schematic representation of the active site of GgNAD-MEβ1. The substitutions in GgNAD-MEβ1 with respect to GgNAD-MEβ2 are highlighted in red.

Recently, we reported that NAD-MEβ1 from the C₄ species G. gynandra and C. angustifolia possess five identically substituted amino acids (V131, S132, H195, V297, and K605) compared with other species and the β2-NAD-MEs of the same species (I131, Y132, Q195, D297, and E605, sequence numbering after alignment to the β1 sequence; Tronconi et al., 2020). To analyze the structural implications of the substitutions in C₄-NAD-MEβ1, we built homology models of the α, β₁, and β₂-isoforms of G. gynandra and A. thaliana. Among the five conserved substitutions in the C₄-NAD-MEβ1 isoforms, only two lead to major changes in the side-chain properties: Y132S and D297V. Both substitutions are located near the catalytic center of the enzyme (Figure 3B). Y132 is directly involved in the catalytic mechanism of the enzyme, donating a proton to the enolpyruvate intermediate to facilitate tautomerization to pyruvate (Chang and Tong, 2003). For that, the hydroxy group needs to be close to the C3 atom of malate. We found that the Y132S substitution almost doubles the distance between the malate C3 and the hydroxy group (3.7–6.8 Å) (Figure 3B). Furthermore, the pK_a of serine is approximately 3 log units larger than that of tyrosine. Thus, serine at position 132 likely cannot donate a proton to enolpyruvate. D297 is not directly involved in the catalytic mechanism but forms a charge-assisted hydrogen bond with K203; the latter is not protonated and accepts the proton from the 2-OH group of malate in the oxidation step (Chang and Tong, 2003). Due to the loss of the polar interaction through the D297V substitution, likely, K203 becomes more mobile, and its pK_a decreases, which weakens its basicity. Hence, we propose that the substitutions Y132S and D297V likely affect the (de-)protonation events that occur during the catalytic cycle of the NAD-ME. This supports the hypothesis that NAD-MEβ1 evolved as a noncatalytic subunit in the C₄ Cleome species G. gynandra and C. angustifolia.

GgNAD-MEβ1 confers regulatory functions on the photosynthetic NAD-ME

The evolutionary changes that occurred in the sequence of NAD-MEβ1 in C₄ Cleome species (Tronconi et al., 2020) let us hypothesize that this non-catalytic subunit was retained because it must have evolved other molecular properties important for its function in C₄ photosynthesis. One such function could be the regulation of the enzymatic activity by metabolic effectors.

To test this hypothesis, we analyzed the influence of selected metabolites on the activity of GgNAD-ME isoforms. From the glycolytic and TCA cycle intermediates analyzed, we found that GgNAD-MEα/β1 activity is strongly enhanced by fructose-1,6-bisphosphate (FBP), PEP, fumarate, and acetyl-CoA (Figure 4A;Supplemental Data Set S1). This agrees with the reported regulation of all A. thaliana NAD-MEs (Tronconi et al., 2010a, 2010b). Also, GgNAD-MEα is strongly activated by fumarate, and GgNAD-MEβ2 is activated by FBP, PEP, and CoA (Figure 4A;Supplemental Data Set S1). As we already showed that FBP, PEP, and CoA bind to the surface of Arabidopsis NAD-MEβ subunit (Tronconi et al., 2010a, 2010b), the activation of GgNAD-MEα/β1 by these effectors provides evidence that GgNAD-MEβ1, although non-catalytic in vitro, is likely a subunit with conserved regulatory functions. Furthermore, GgNAD-MEα/β1 shows specific inhibition by 3-phosphoglycerate (3PG) and activation by oxaloacetate (OAA) (Figure 4A;Supplementary Data Set S1).

Figure 4

Regulatory properties of recombinant GgNAD-MEs. A, Modulation of the activity of GgNAD-ME entities by selected metabolites. The values represent the mean ± sd (standard deviation); n = three independent replicates. Enzymatic activities were measured in the absence or presence of each effector. Malate concentration was kept at a value of K_m/3 for each enzyme. The results are presented as the log₂ FC (FC, fold change = enzymatic activity in the presence of each effectors/enzymatic activity measured in the absence of the effectors). FC values <0.7 (log₂ < −0.5) or >1.3 (log₂ > 0.4) were significantly different to 1 (P <0.05, one-way ANOVA with Holm–Sidak method) and considered as an inhibition (left bars) or activation (right bars) of the enzymatic activity, respectively. Red arrows indicate metabolites that differentially regulated GgNAD-ME α/β1 activity (see Supplemental Data Set S1 for the complete data set). B, Cartesian plane representation of the kinetic parameters of recombinant GgNAD-MEα/β1 (Ggα/β1), GgNAD-MEα/β2 (Ggα/β2), ThNAD-MEα/β1 (Thα/β1), ThNAD-MEα/β1 (Thα/β2), AtNAD-MEα/β (Atα/β), and the interspecies chimera AtNAD-MEα/GgNAD-MEβ1 (Atα/Ggβ1) in the absence and presence of aspartate (see Supplemental Table S3 for the complete data set). Mean values of K_m malate and k_cat are shown as an ordered pair (K_m, k_cat) with error bars indicating ± se. The blue arrows connect the kinetic parameter values in the absence of aspartate (gray squares) and those obtained in the presence of 5-mM aspartate (white squares). The red arrow indicates the replacement of the β-subunit, which is associated with the Atα-subunit in the interspecies chimera Atα/Ggβ1. The gray area separates enzymes with k_cat/K_m values <100 s⁻¹*mM⁻¹. All points on the boundary line have a constant k_cat/K_m value (= 100). The boundary line was generated from the M–M equation in the form rate =  (k_cat/K_m)E₀[S]/(1+[S]/K_m) at E₀ = 1 (arbitrary units) and [S]  = 0.5  mM (S: malate), setting rate as dependent and K_m as independent variable, respectively, varying K_m from 0.02 to 8 mM.

From the amino acids tested, glycine and aspartate specifically stimulated GgNAD-MEα/β1 activity with apparent constants of activation (A₅₀) of 6.2 mM for glycine and 3.8 mM for aspartate. Aspartate has a stronger inhibitory effect on the activity of GgNAD-MEβ2 and GgNAD-MEα/β2 than glycine and does not modify GgNAD-MEα activity (Figure 4A;Supplemental Data Set S1). Interestingly, in GgNAD-MEβ2, aspartate behaved as a competitive inhibitor of malate (Supplemental Data Set S1; Supplemental Table S3), suggesting that it binds to the active site of the β-subunit.

Given that glycine and aspartate were not previously described as regulators of any plant NAD-ME, we compared their effects on the activities of GgNAD-MEα/β1, AtNAD-MEα/β, and the interspecies chimera AtNAD-MEα/GgNAD-MEβ1. Glycine had no effect on the activities of AtNAD-MEα/β and AtNAD-MEα/GgNAD-MEβ1 (Supplemental Data Set S1). However, aspartate also activates AtNAD-MEα/GgNAD-MEβ1 (Figure 4B;Supplemental Table S3) with an A₅₀ value of 4.5 mM. In the presence of 5-mM aspartate, the catalytic efficiencies of GgNAD-MEα/β1 and AtNAD-MEα/GgNAD-MEβ1 increase up to three-fold due to increases in both malate affinity and turnover rate (Figure 4B;Supplemental Table S4). Moreover, in the absence of aspartate, the affinity for malate of AtNAD-MEα/GgNAD-MEβ1 is three-fold higher than that of AtNAD-MEα/β (Figure 4B;Supplemental Table S3). Aspartate decreases the enzymatic efficiency of AtNAD-MEα/β by a similar magnitude as in GgNAD-MEβ2 and GgNAD-MEα/β2 (Figure 4B;Supplemental Table S3).

All these results indicate that (1) the GgNAD-Meβ1 subunit is involved in the specific kinetic properties of GgNAD-Meα/β1, high affinity for malate and turnover rate, and activation by aspartate; (2) in the β-subunits, aspartate most probably binds to the active site and exerts different regulatory effects depending on the evolutionary adaptations of the subunits; (3) despite the evolutionary divergence between species, the interaction of the NAD-MEα- and β-subunits is largely conserved; and (4) the association of the noncatalytic GgNAD-MEβ1 subunit to the paralogous α-partner of Arabidopsis renders a NAD-ME entity with similar properties as GgNAD-MEα/β1.

GgNAD-MEβ1 stabilizes an associated NAD-MEα subunit

To scrutinize how the noncatalytic GgNAD-MEβ1 subunit influences the activity of an associated GgNAD-MEα subunit, we performed molecular dynamic simulations of four dimeric NAD-ME with aggregate simulation times of 2.5 μs each (Supplemental Figure S5). We combined a “catalytic” NAD-MEα subunit with “effector” subunits α, β1, and β2, resulting in the dimers GgNAD-MEα/α, GgNAD-MEα/β1, GgNAD-MEα/β2, and ThNAD-MEα/β1. Then, we subjected the conformational ensembles to a constraint network analysis (CNA) to reveal the hierarchy of structural rigidity of the dimers. We focused on the effects of the “effector” subunits on the “catalytic” NAD-MEα subunit.

A comparison of the chemical potential energies due to noncovalent bonding per residue (see Eq. 1, “Material and Methods”) for the “catalytic” NAD-MEα subunits showed that the GgNAD-MEβ1 subunit has a more stabilizing influence on the associated GgNAD-MEα subunit than any of the other “effector” subunits (Figure 5, A and B; Supplemental Figure S6). We found that the more stabilizing influence of GgNAD-MEβ1 originates from forming more rigid contacts (rc) across the interface to the GgNAD-MEα subunit than any of the other “effector” subunits (Figure 5C;Supplementary Figure S6B). This stabilizing effect spreads from the interface region to the active site (Figure 5, B and C; Supplemental Figure S6A). The higher structural stability of the active site likely improves the binding of the substrate and cofactor or their orientation, leading to the higher catalytic efficiency of GgNAD-MEα/β1 with respect to all other NAD-ME entities analyzed (Figure 2C). The same mechanism may lead to the increased catalytic efficiency of the interspecies chimera AtNAD-MEα/GgNAD-MEβ1 with respect to the AtNAD-MEα/β dimer (Figure 4B;Supplemental Table S3).

Figure 5

Molecular dynamic simulations of the NAD-MEα subunit with “effector” subunits. A, The difference profile of the chemical potential energy due to noncovalent bonding per residue of the GgNAD-MEα subunits of a GgNAD-MEα/α homodimer and a GgNAD-MEα/β1 heterodimer. Values <0 indicate a stabilizing effect of the β1-subunit on the α-subunit with respect to the homodimer, values >0 indicate a destabilizing effect. Above the plot, a schematic representation of the enzymes’ secondary structure is shown, with markers below highlighting residues in the subunit interface and markers above highlighting residues in the active site. B, The difference profile from A mapped onto the “catalytic” α-subunit of a GgNAD-MEα/β1 heterodimer (with bound NAD⁺ and L-malate), with the “effector” subunit being depicted as a white ribbon. Blue (red) indicates more (less) structurally stable regions than in the GgNAD-MEα/α homodimer. Corresponding representations of all analyzed systems are shown in Supplemental Figure S6A. C, Depiction of rc strength between amino acid pairs in the GgNAD-MEα/α homodimer (left) and the GgNAD-MEα/β1 heterodimer (right) with bound NAD⁺ and L-malate. A rc denotes that two residues are part of a structurally rigid cluster in the structure. Blue (yellow) indicates more (less) prevalent rc. The “catalytic” α subunit is depicted as helices and sheets at the top, and the “effector” subunit as ribbons below. In the GgNAD-MEα/β1 heterodimer, there are more rc across the interface than in the GgNAD-MEα/α homodimer.

Furthermore, we investigated how specific substitutions in the NAD-MEβ1 subunit in C₄ Cleome may stabilize an associated NAD-MEα subunit. For this, we applied an ensemble-based perturbation approach implemented in constrained CNA, perturbing selected residues in GgNAD-MEβ1 versus the corresponding residues in GgNAD-MEβ2. We applied the approach on four substitution sites that we previously showed to be identically substituted in C. angustifolia and G. gynandra NAD-MEβ1 with respect to other NAD-MEβ enzymes: I131V, Y132S, Q195H, and D297V (Tronconi et al., 2020). We found that, although the I131V substitution shows no impact on the rigidity of GgNAD-MEβ1, all other substitutions reduce the stability around the inactive catalytic site (Supplemental Figure S7). This destabilization, however, remains localized in the NAD-MEβ1 subunit. Interestingly, the Q195H substitution additionally exhibits a stabilizing effect in the interface region, which also spreads to the bound NAD-MEα subunit (Supplemental Figure S7C). These results suggest that H195 in C. angustifolia and G. gynandra NAD-MEβ1 confer stabilization of the associated subunit.

Immunolocalization studies

Our results indicate that in C₄ Cleome, a housekeeping NAD-MEα/β2 isoform and a photosynthetic NAD-MEα/β1 (C₄-NAD-ME) isoform coexist in photosynthetic tissues and are formed by a differential combination of subunits. This gives rise to two hypotheses regarding the cellular location of both enzymatic entities in photosynthetic tissues of C₄ Cleome. Either housekeeping and C₄-NAD-ME are separated by cell type, or the housekeeping NAD-ME is found in every cell type, and only C₄-NAD-ME is exclusively located in BSCs. To answer this question, we conducted an immunolocalization analysis using leaves of G. gynandra and specific antibodies raised against each of the three NAD-ME subunits. We found high-density labeling of GgNAD-MEα and GgNAD-MEβ2 subunits in mitochondria of mesophyll cells (MC) (Figure 6). The use of antibodies against GgNAD-MEβ1 rendered few spots, indicating either very low occurrence of the NAD-MEβ1 subunit in MC mitochondria or just background labeling. Mitochondria of BSC presented immunolabeling of all three subunits and at much higher density than those in MC (Figure 6;Supplemental Figure S8). These results indicate that the housekeeping and photosynthetic NAD-ME entities coexist in BSC mitochondria of the C₄ Cleome species (Figure 7).

Figure 6

Immunolocalization analysis in G. gynandra leaves. Transmission electron micrographs of bundle sheath (A–C) and mesophyll (D–F) cells of G. gynandra showing mitochondria with immunogold labeling (black particles) of GgNAD-MEα (A and D), GgNADMEβ1 (B and E), and GgNAD-MEβ2 (C and F). *mitochondrion. Bars, 500 nm (A–C); 100 nm (D–F).

Figure 7

Schematic representation of the localization of C₄ Cleome G. gynandra NAD-ME entities in different types of cells. GgNAD-MEα/β2 is present in mitochondria of every cell besides BSCs where it is involved in malate respiration as an associated enzyme of the TCA cycle. GgNAD-MEα/β1 is present in mitochondria of BSC where it performs the C₄ decarboxylase activity delivering CO₂ for the Calvin Benson cycle.

Discussion

A straightforward rationale is to assume that if a gene was recruited into the C₄ pathway, it would be expressed at higher levels than in C₃ species and its transcript would mostly accumulate during the day (Hibberd and Covshoff, 2010; Brautigam et al., 2011; Pick et al., 2011; Christin et al., 2013a, 2013b; Moreno-Villena et al., 2018). This is the case for most of the C₄ genes (Brautigam et al., 2011; Pick et al., 2011). Nevertheless, our comparative quantitative transcriptional analysis revealed that all NAD-ME genes are more highly expressed in C₄ Cleome leaves than their homologs in C₃ Cleome leaves. Therefore, it is not possible to unambiguously identify a concrete C₄-NAD-ME candidate gene by the expressional analysis, which in part explains why, until now; a specific C₄-NAD-ME has not been identified at the molecular level.

Here, using a combination of biochemical and proteomic analyses, we identify different NAD-ME isoforms in C₃ and C₄ Cleome organs, which arise through a differential combination of subunits previously shown to be selectively adapted through C₄ evolution (Tronconi et al., 2020). We found that non-photosynthetic and photosynthetic organs of C₃ Cleome T. hassleriana possesses two NAD-ME isoforms, ThNAD-MEα/β1 and ThNAD-MEα/β2, with kinetic properties similar to those of A. thaliana NAD-MEs. These results indicate that ThNAD-MEα/β1 and ThNAD-MEα/β2 represent housekeeping enzymes that participate in mitochondrial malate respiration. The paralogs β1 and β2 genes were possibly fixed by genetic drift (Tronconi et al., 2020), most probably resulting in redundancy of function for ThNAD-MEα/β1 and ThNAD-MEα/β2 in the C₃ species. In the C₄ Cleome species, NAD-MEα/β2 is responsible for the housekeeping function.

We found that NAD-MEα/β1 is present exclusively in leaves of the C₄ species G. gynandra and C. angustifolia. Here, NAD-MEα/β1 exists in BSC mitochondria together with the NAD-MEα/β2 housekeeping heteromer. GgNAD-MEα/β1 has a high affinity for malate and high catalytic efficiency, and its activity is enhanced in the presence of aspartate, an intermediate of the NAD-ME subtype C₄ pathway. Aspartate is formed in the mesophyll cells from OAA, the product of the first carboxylation step, and is used to shuttle the prefixed CO₂ to BSC mitochondria (Ludwig, 2016a, 2016b). Aspartate moves from mesophyll cells to the BSC by diffusion; in this case, aspartate would be expected to be maintained at quite low levels in the BSC cytosol. Most probably, active uptake of aspartate into the mitochondria would aid the diffusion of aspartate from mesophyll to BSC. In fact, the existence of uncoupling proteins 1 and 2 (UCP1 and UCP2), two mitochondrial protein carriers that catalyze a highly efficient electroneutral exchange of aspartate and glutamate in Arabidopsis (Monne et al., 2018), supports this idea. In BSC mitochondria, aspartate is then converted to malate. In G. gynandra, the aspartate aminotransferase recruited into the C₄ pathway is located in the mitochondria (Sommer et al., 2012). Apart from that, mitochondrial aspartate can be used to speed up the C₄ cycle through activation of C₄-NAD-ME. Due to the mentioned kinetic and regulatory properties of GgNAD-MEα/β1 and the fact that this enzymatic entity is only found in leaves of C₄ Cleome, we conclude that GgNAD-MEα/β1 is the C₄-NAD-ME photosynthetic isoform.

Our extensive enzymatic analysis indicated that NAD-MEβ1 of both C₄ species, G. gynandra and C. angustifolia, are not active. Comparison of closely related β-NAD-ME sequences indicated that NAD-MEβ1 of the C₄ species accumulated many amino acid changes during evolution (Tronconi et al., 2020). We found that the changes of Y132 and D297, which participate in the catalytic mechanism, are most probably responsible for the catalytic inactivity of Ca/GgNAD-MEβ1. The loss of function of the active site could have led to the generation of an allosteric site for aspartate, a C₄ compound structurally similar to malate. Our results indicate that aspartate competitively inhibits NAD-MEβ2 (Supplemental Tables S3 and S4), suggesting that it binds effectively to the active site. Liu et al. (2002) showed that NAD-ME from Ascaris suum catalyzes the oxidative decarboxylation of aspartate to pyruvate and ammonia through a mechanism that likely mimics that of malate. Thus, the evolution from an active to an allosteric site would have required few amino acid substitutions.

We also found that the Q195H substitution in the NAD-MEβ1 has a stabilizing effect in the dimer interface region (Supplementary Figure S7C). The stabilizing influence of NAD-MEβ1 on the associated NAD-MEα most probably enables the differential assembly of the subunits to form independent protein entities acting in different metabolic pathways in the same cell compartment. The flux of malate through these metabolic pathways is likely regulated by the action of aspartate as we found that this C₄ intermediate highly enhances the efficiency of C₄-NAD-ME but at the same time decreases the efficiency of the respiratory isoform (Supplemental Data Set S1 and Supplemental Table S3). The contrasting behavior of the NAD-MEβ subunits in C₄ Cleome points to a specialized role of NAD-MEβ1 in the C₄ isoform.

Phylogenetic analysis and differential types of Kranz anatomy in G. gynandra and C. angustifolia indicated that they have two different C₄ origins (Feodorova et al., 2010; Koteyeva et al., 2011; Bayat et al., 2018). Although these C₄ species display differences in leaf architecture and physiology (Koteyeva et al., 2011), the recruitment of NAD-ME genes into the C₄ pathway followed a similar molecular mechanism. After fixing the pair of paralog β genes, one copy was co-opted to play a role in C₄ photosynthesis through parallel amino acid substitutions. We postulate that the noncatalytic subunit in Cleome C₄ species, NAD-MEβ1, was evolutionarily fixed because it enabled the formation of a NAD-MEα/β1 heteromer with suitable biochemical properties to fulfill a role in the C₄ photosynthetic pathway, such as high catalytic efficiency and activation by aspartate.

In summary, we report how the neofunctionalization of a heteromeric enzyme can occur through the duplication and subsequent adaptation of one subunit, while the other subunit is shared between the original and new enzyme. Our work paves the way to elucidate the molecular mechanism underlying protein recruitment for novel functions, especially for proteins with heteromeric composition. Specifically, our work shows that a C₄-exclusive NAD-ME isoform exists in plant cell mitochondria, which is formed through the differential assembly of protein subunits that accumulated adaptive mutations during evolution. This discovery will now aid in introducing C₄ traits into C₃ plants (von Caemmerer and Furbank, 2016; Ermakova et al., 2020; Lin et al., 2020). Ongoing projects can now rationally design the incorporation of a NAD-ME with C₄ characteristics by introducing the NAD-MEβ1 of C₄ Cleome or modifying a preexisting extra NAD-MEβ copy.

Materials and methods

Plant growth conditions

Gynandropsis gynandra and T. hassleriana were grown under greenhouse conditions in commercially available soil mixed with Cocopor (Stender, Schermbeck, Germany). G. gynandra was germinated in the dark for 3–5 days, while T. hassleriana needed ∼14 days for germination under a long day (16 h/8 h) light regime. Natural light was supplemented with regular filament lamps mounted 1.2 m above the surface, ensuring a minimum of 80 µmol m⁻² s⁻¹ PAR. Seedlings were transferred from germination pots to single clay pots after full formation of the cotyledons. Plants were watered from above as needed and single pot soil was premixed with 3 g L⁻¹ Osmocote (Scotts Deutschland GmbH, Nordhorn, Germany) as fertilizer. Arabidopsis thaliana (ecotype Columbia-0) was grown in soil as described before (Hüdig et al., 2015). For immunolocalization studies, plants were grown in 20 L pots, watered daily to avoid drought stress, and fertilized weekly as described (Sage et al., 2013). Plants were exposed to direct sunlight most days, giving light intensities of 1,600 μmol m⁻² s ⁻¹.

Transcriptional analysis of NAD-ME genes

Total RNA from samples were isolated using the RNeasy Mini Kit (Qiagen, Hilden, Germany) from pooled (n = 5) plant leaf material of G. gynandra and T. hassleriana. Genomic DNA was removed by DNAse treatment (Ambion Inc., Austin, TX, USA), according to the manufacturer’s instructions. cDNA synthesis was performed using RevertAid H Minus Reverse Transcriptase (Thermo Fisher Scientific, Darmstadt, Germany) with oligodT primers according to the manufacturer’s instructions. Primers for quantitative real-time polymerase chain reaction (qRT-PCR; Supplemental Table S5) were designed to amplify a PCR product of 130–210 bp length in the region 450–750 bp from the 3′-end of the coding sequence of the mRNA transcript. The respective Actin2 homolog was chosen from the nearest orthologs to A. thaliana Actin2. Primer efficiency was determined with regards to primer concentration and annealing temperature, and was used to normalize C_t values. Primer specificity was ensured using cloned coding sequences from plasmids as test templates and agarose gel visualization of the respective final RT-qPCR products. RT-qPCR was run on StepOne Plus (Thermo Fisher Scientific, Darmstadt, Germany) equipped with StepOne Software version 2.2.2 using KAPA SYBR FAST qPCR Kit Master Mix (2X) ABI Prism (Kapabiosystems, Boston, MA, USA) in a 15 µL reaction. RT-qPCR amplification was performed at: 95°C for 3 min, 40 cycles at 95°C for 3 s and 65°C for 20 s, followed by a standard melting curve.

Cloning

Coding sequences of the NAD-ME subunits of G. gynandra and T. hassleriana were amplified from leaf cDNA using Phusion Polymerase (Thermo Fisher Scientific, Darmstadt, Germany) and specific primers (Supplemental Table S4). The amplified fragments were sub-cloned into pCR-TOPO-BluntII (Thermo Fisher Scientific, Darmstadt, Germany) and sequenced using EZ-seq sequencing services (Marcorgen Europe, Amsterdam, The Netherlands). The generated TOPO vectors were used as source for the generation of cohesive-end insert fragments for restriction site cloning or as templates for further PCR-fragment constructions according to Gibson (2011). The vector pET16b (Merck, Darmstadt, Germany) was used for heterologous expression in E. coli. The pET16 constructions produce NAD-MEs containing a His•Tag (10 histidine residues) sequence fused to the N-terminus for its purification through affinity chromatography. When using classic restriction enzyme cloning, the pET16b vector was linearized with restriction enzymes (NdeI/XhoI/BamHI) as well as the respective PCR products to generate the cohesive ends, and circular plasmids were reconstituted with T4 Ligase (Thermo Fisher Scientific, Darmstadt, Germany). If Gibson assembly was used, the primers were designed with a 20 bp overlap with the destination vector. Purified Gibson-assembly PCR-fragments of NAD-ME coding sequences were cloned into pET16b linearized with BamHI. Gibson isothermal assembly was performed according to Gibson (Gibson, 2011) with reagents received from Thermo Fisher Scientific (Darmstadt, Germany) and New England Biolabs (Ipswich, MA, USA). Final plasmids were sequenced and amplified using E. coli DH5α cells (Thermo Fisher Scientific, Darmstadt, Germany).

For the co-expression of NAD-ME subunits, the cDNA fragments corresponding to the mature NAD-MEα subunits of G. gynandra and T. hassleriana were cloned in the pET32 vector (Novagen) using specific restriction sites. The restriction sites NcoI/SalI were used for GgNAD-MEα and BamHI/SalI for GgNAD-MEβ1 and -β2. The pET32-NAD-MEα constructions produce NAD-MEα containing Trx•Tag (109 aa of thioredoxin) and His•Tag sequences fused to the N-terminus of NAD-ME. On the other hand, the cDNA fragments coding for the mature NAD-MEβ subunits were cloned in the pET29 vector (Novagen, Vadodara, India), which adds a C-terminal His•Tag sequence. The pET29-NAD-MEβ constructions express the mature NAD-MEβs without the fusion vector-coded sequence because the cloned cDNA has a stop codon in its 5′end. For the co-expression of GgNAD-MEβ1 and -β2, alternatively one β-subunit-coding sequence was cloned in the pET32 vector.

The coding sequence of the mature C. angustifolia NAD-MEβ1 subunit was synthetized using the BioCat commercial gene synthesis service (Heidelberg, Germany) and cloned into pET32b via NcoI/BamHI restriction sites. The resulting plasmid was sequenced and verified using the Marcogen sequencing service.

Heterologous expression and purification of recombinant NAD-MEs

Recombinant single NAD-ME subunits were produced with the T7-polyperase IPTG-inducible expression system (Thermo Fisher Scientific, Darmstadt, Germany). For the production of His-tagged GgNAD-MEα and ThNAD-Meα, we used E. coli BL21 (DE3) cells (Merck, Darmstadt, Germany) (from now on denoted as protocol 1, P1) and for production of His-tagged GgNAD-MEβ1, GgNAD-MEβ2, CaNAD-MEβ1, ThNAD-MEβ1, and ThNAD-MEβ2, we used E. coli ArcticExpress (DE3) cells (Agilent Technologies, Santa Clara, USA) (from now on denoted as protocol 2, P2). Chemically competent E. coli cells were transformed with 50 ng of the expression vector and grown on LB agar plates with antibiotics according to the manufacturer’s instructions.

For the co-expression of NAD-ME subunits, E. coli BL21 (DE3) cells were simultaneously transformed with the pET29 and pET32 constructions. The cells were selected on LB-agar plates supplemented with 100 µg/mL ampicillin (pET32 selection agent) and 50 µg/mL kanamycin (pET29 selection agent). For protein expression, 800 mL liquid cell cultures containing 100 µg/mL ampicillin and 50 µg/mL kanamycin were inoculated with freshly grown overnight culture for 3 h to an OD₆₀₀ of 0.6. Protein expression was induced with 1 mM IPTG. Cells were harvested after 16 h at 16°C by centrifugation for 10 min at 4°C and 6,000 g, and stored at −20°C until use. The interspecies chimera AtNAD-MEα/GgNAD-MEβ1 was obtained by cotransformation of E. coli BL21 (DE3) with pET32-AtNAD-MEα (Tronconi et al., 2008) and pET29-GgNAD-MEβ1.

Induction of protein expression and purification of recombinant NAD-MEs

For induction of protein expression, 400 mL liquid cell cultures containing 100-µg/mL ampicillin (P1) or 100-µg/mL ampicillin and 20-µg/mL gentamicin (P2) were inoculated with freshly grown over-night culture, grown at 37°C for 2.5 h to an OD₆₀₀ of 0.5 (P1) or at 37°C for 3 h to an OD₆₀₀ of 0.8 (P2). Protein expression was induced with 1-mM IPTG. Cells were harvested after 20 h growth in the dark at 37°C (P1) or 3 days at 12°C (P2) by centrifugation for 10 min at 4°C and 6,000 g, and stored at −20°C until use for a maximum of 3 months. His-tagged single proteins were purified using the immobilized metal ion affinity chromatography principle. Escherichia coli cells were resuspended in freshly prepared ice-cold lysis buffer (0.2 M NaCl, 20 mM Tris–HCl pH 8.0, 5 mM imidazole, 2 mM phenylmethylsulfonyl fluoride (PMSF), 1 spatula tip lysozyme), incubated for 10 min on ice, and sonified for 4 min in 30 s intervals. After centrifugation (15 min, 15,000g, 4°C), supernatant was loaded on a 900 µL nickel-nitrilotriacetic acid (Ni-NTA) agarose (Qiagen, Hilden, Germany) column. Three washing steps with incrementing imidazole concentration were performed (0.2 M NaCl, 20 mM Tris–HCl pH 8.0, and 5, 40, 80 mM (P1) or 5, 40, 60 mM (P2) imidazole) at 4°C with 10 column volumes each until the purity of the respective protein was achieved. Proteins were eluted in 2 mL elution buffer (0.2 M NaCl, 20 mM Tris–HCl pH 8.0, 0.5 M imidazole). The elution of NAD-MEα, -β1, and -β2 proteins was verified by Coomassie-stained sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS–PAGE). The expected molecular masses were ∼67.5, 66.3, and 66.4 kDa for NAD-MEα, -β1, and -β2, respectively.

For purification of co-expressed proteins, cells were resuspended in 80 mL buffer A (20 mM Tris–HCl, pH 7.9, 5 mM imidazole, and 2 mM PMSF), sonicated, and centrifuged for 10 min at 7,000 g at 4°C. The supernatant was loaded onto a Ni-NTA agarose (Qiagen, Hilden, Germany) column previously equilibrated with buffer A. The co-expressed proteins were eluted with buffer A containing 200-mM imidazole. The co-elution of both NAD-MEα and -β proteins was verified by Coomassie-stained SDS–PAGE. The expected molecular masses were 75–82 kDa for the pET32-expressed proteins (58–65 kDa corresponding to the NAD-ME + 17 kDa the His-Trx•Tag) and 58–65 kDa in the case of NAD-MEs expressed in the pET29 vector. The equimolar ratio of the co-purified proteins was estimated by densitometric analysis of the bands obtained after SDS–PAGE of the eluted fractions.

The co-purified fusion NAD-MEs were incubated with thrombin protease (1:100) for 2 h at 15°C to remove the N-terminus encoded by the pET32 vector. The proteins were further purified using a gel chromatography Sephadex G-50 column equilibrated with buffer B (50 mM MES pH 6.5, 5 mM MnCl₂, 5 mM dithiothreitol (DTT), and 20% (v/v) glycerol). Purified enzymes were stored at −80°C in buffer B with 50% glycerol, as previously described (Tronconi et al., 2008).

Protein quantification and gel electrophoresis procedures

SDS–PAGE was used to monitor protein purification according to Laemmli (1970). Gradient native PAGE (Novex 4%–12% Tris–Glycine gels, Invitrogen, Thermo Fisher Scientific) was used to separate native protein complexes from plant tissue by omitting β-mercaptoethanol and SDS from any buffer used for sample loading, pouring, or running PAGE gels. Gels were run at 4°C and 100 V. Proteins in gels were visualized using Coomassie brilliant blue staining. Protein concentration was determined using a Bio-Rad Protein Assay (Bio-Rad, Hercules, USA) or, alternatively, a simplified amido black 10B precipitation method (Schaffner and Weissmann, 1973) using total serum protein as standard.

Preparation of protein extracts for native PAGE and in-gel NAD-ME activity assays

Plant tissue was freshly harvested from greenhouse grown G. gynandra, T. hassleriana, C. angustifolia, and A. thaliana and ground to fine powder in liquid nitrogen. Soluble proteins were extracted by adding 3 µL per milligram fresh weight ice-cold extraction buffer (0.1 M Tris–HCl pH 8.0, 0.1 M NaCl, 0.5% (v/v) Triton X-100, 2 mM PMSF, 10 mM DTT, 0.1% (w/v) polyvinylpyrrolidone 40), followed by vigorous vortexing and incubation for 20 min on ice. Soluble protein was separated by centrifugation (10,000g, 10 min, 4°C) and supernatant was transferred to a new tube. Proteins were separated in native PAGE run at 100 V at 4°C.

After electrophoresis, gels were assayed for NAD(P)-ME activity. For NAD-ME and NADP-ME activity assay, gels were first incubated in 50 mM MOPS-NaOH (pH 6.8) at room temperature for 30 min. Rebuffered native PAGE gels were incubated in 50 mM MOPS-NaOH (pH 6.8), 0.05% (w/v) nitro blue tetrazolium, 150 μM PMS, 10 mM MnCl₂, 5 mM NAD, or 5 mM NADP, and either 20 mM, 100 mM or no malate. Gels were documented after initial clearance in ddH₂O and gel slices with NAD-ME activity were processed for protein identification by mass spectrometry.

Mass spectrometry

Slices of PAGE samples for mass spectrometry were essentially processed as described previously (Brenig et al., 2020). Briefly, protein-containing bands were washed, reduced by DTT, alkylated with iodoacetamide, and digested with trypsin. Resulting peptides were extracted from the gel and cleaned up using solid-phase extraction (HLB µElution plate, Waters) and finally resuspended in 0.1% trifluoroacetic acid and analyzed by liquid chromatography coupled mass spectrometry. Peptide separation was carried out on an Ultimate 3000 rapid separation system (Thermo Fisher Scientific) on 25-cm long C18 columns using a 54-min separation gradient (2 h for G. gynandra leaf band migration level 1 sample 2) essentially as described (Poschmann et al., 2014). Subsequently, separated peptides were injected via a nano-electrospray interface into the mass spectrometer operated in data-dependent positive mode and searches were carried out in the proteome discoverer framework as follows. Gynandropsis gynandra and T. hassleriana samples were analyzed on a QExactive Plus mass spectrometer. Briefly, full scans were carried out at a resolution of 70,000 (scan range 200–2,000 m/z, profile mode, automatic gain control target value 3E6, maximum injection time 50 ms). Subsequently, up to 20 two- to five-fold precursors were selected (4 m/z isolation window), fragmented by higher-energy collisional dissociation (HCD) and analyzed at a resolution of 17,500 (available scan range 200–2,000 m/z, centroid mode, automatic gain control target value 1E5, maximum injection time 50 ms), the dynamic exclusion was 10 s. Parameters for the G. gynandra leaf band migration level 1 sample 2 were as indicated above with the following modifications: MS1 scan range, 350–2,000 m/z; MS1 maximum injection time, 80 ms; isolation window, 2 m/z; MS2 maximum injection time, 60 ms; top 10 method; 100 s dynamic exclusion.

The C. angustifolia proteins were analyzed on a Lumos Fusion hybrid instrument (Thermo Fisher Scientific) with the following parameters: MS1 (Orbitrap) resolution 120,000, MS1 scan range 200–2000 m/z, MS1 automatic gain control target 400,000, MS1 maximum injection time 60 ms, 1.6 m/z isolation window (charge states 2–7 included), HCD fragmentation, MS2 (linear ion trap) scan rate rapid, MS2 automatic gain control target 400,000, MS2 maximum injection time 150 ms, cycle time 2 s, 60 s dynamic exclusion.

Peptide and protein searches were carried out within the Proteome Discoverer 2.4.1.15 framework (Thermo Fisher Scientific). MS Amanda 2.0 was used as a search engine followed by a fixed value peptide spectrum matches (PSMs) validator node and peptides and proteins accepted at a false discovery rate of 1%. PSMs were reported only for unambiguously assignable peptides and if a minimum of two PSMs were identified for one protein. For MS Amanda, the following settings were applied: tryptic cleavage specificity, maximum of two missed cleavages, carbamidomethylation at cysteines as fixed, and protein N-terminal acetylation and methionine oxidation as variable modifications, MS1 tolerance 5 ppm, MS2 tolerance 0.02 Da (0.4 Da for C. angustifolia proteins). Acquired spectra were matched against transcriptome-derived G. gynandra (32,832 entries) and T. hassleriana (28,718 entries) protein databases based on the transcriptome datasets of Kulahoglu et al. (2014). For the identification of C. angustifolia proteins, a custom database was created using the G. gynandra transcriptome and adding the known NAD-ME coding sequences.

Kinetic analysis

Purified recombinant NAD-MEs were assayed for the oxidative decarboxylation of malate using a standard reaction mixture containing 50 mM MES-NaOH pH 6.2–6.4, 10 mM MnCl₂, 4 mM NAD, 20-mM malate in a final volume of 0.5 mL. The reaction was started by the addition of malate. An extinction coefficient (ε) of 6.22 mM⁻¹ cm⁻¹ for NADH at 340 nm was used in the calculations. Initial-velocity studies were performed by varying the concentration of one of the substrates around its K_m value while keeping the other substrate concentrations at subsaturating or saturating levels. All kinetic parameters were calculated at least by triplicate determinations and adjusted to non-lineal regression using free concentrations of all substrates. One unit (U) is defined as the amount of enzyme that catalyzes the formation of 1 mmol of NADH per min under specified conditions.

When testing the ability of metabolic intermediates as possible modulators (inhibitors or activators) of the enzymatic activity, NAD-ME activity (1–3  µg of protein) was measured in the absence or presence of 0.5 or 5.0 mM of this metabolite, at a saturating NAD concentration (4 mM) and sub-saturating levels of malate (1 mM for the housekeeping isoforms and 0.1 mM for GgNAD-MEα/β1), which were kept at one-third of the K_m value for each enzyme (K_m/3). The apparent activation constant (A₅₀) values were obtained by varying the concentration of activator between 0.1 and 30 mM, while keeping the NAD concentrations at a saturating level (4 mM) and malate at nonsaturating concentrations (∼3 mM for the housekeeping isoforms and ∼0.3 mM for GgNAD-MEα/β1).

Gel filtration chromatography

Molecular masses of recombinant GgNAD-MEs were evaluated by gel-filtration chromatography on an ÄKTA purifier system (GE Healthcare) using a Biosep-SEC-S3000 column (Phenomenex). The column was equilibrated with 50 mM MES-NaOH, pH 6.5, and calibrated using molecular mass standards and calibrated using molecular mass standards (Sigma-Aldrich, Merck, Darmstadt, Germany). The samples and standards were applied in a final volume of 100 µL at a constant flow rate of 0.5 mL/min. All the experiments were performed in triplicate.

Generation of comparative models

We used SWISS-MODEL (Waterhouse et al., 2018) to build comparative models of the α, β1, and β2 isoforms of GgNAD-ME as well as the α and β₁ isoform of ThNAD-ME. All models were built using PDB ID 1PJ3 (Tao et al., 2003) as a template structure (sequence identities: 42.0%–42.8%, coverage: >88%). After superimposing the protein structures, the coordinates of NAD⁺, L-malate, and Mg²⁺ in PDB ID 1PJ2 (Tao et al., 2003) were copied to the comparative models. Crystal waters within 4.0 Å of any of the placed molecules were taken from PDB ID 1PJ2, too. Molecules bound to the template outside the active site were not considered. Using the protein preparation wizard of MAESTRO (Schrödinger Release 2018-1: Maestro; Schrödinger, LLC, New York, NY, 2018), we capped the termini and protonated the structures according to pH 7.0 using Epik (Schrödinger Release 2018-1: Epik; Schrödinger, LLC, New York, NY, 2018). Unusual protonation states (protonated His or protonated Asp and Glu) were visually inspected; the protons were afterward removed. The lysine residue interacting with the 2-OH group of L-malate was set to a deprotonated, uncharged state to simulate the state of the enzyme before the oxidation of L-malate. All hydrogen were minimized using the OPLS3e force field while keeping all heavy atoms fixed. Afterward, crystal water molecules from PDB ID 1PJ3 further away than 2.2 Å away from L-malate, NAD⁺, and Mg²⁺ were placed in the structures. The protonated structures of the monomeric isoforms were used to build the GgNAD-MEα/α homodimer as well as the GgNAD-MEα/β1, GgNAD-MEα/β2, and ThNAD-MEα/β1 heterodimers by superimposing the monomers onto 1PJ3.

Molecular dynamic simulations

tLEaP from the Amber20 software package (AMBER 2020; University of California, San Francisco, 2020) was used to prepare the complexes for the molecular dynamics (MD) simulations. The complexes were solvated in an octahedral box of TIP3P water (Jorgensen et al., 1983), leaving at least 11.0 Å to any edge of the box. Random water molecules were replaced by sodium to neutralize the charges. ff14SB (Maier et al., 2015) was used for the proteins. For L-malate and NAD⁺, GAFF2 parameters were used, and atomic point charges were derived using the restrained electrostatic potential (RESP) procedure (Bayly et al., 1993) from structures optimized with Gaussian 16 (Gaussian 16 Rev. A.03; Gaussian, Inc., Wallingford CT, 2016) at the HF/6-31G* level of theory. The Mg²⁺ cations in the catalytic site were modeled employing a cationic dummy atom model (Jiang et al., 2015).

All MD simulations were performed using the graphics processing unit version of pmemd (Le Grand et al., 2013; Salomon-Ferrer et al., 2013). The Langevin thermostat (Pastor et al., 1988) was used for temperature control with a collision frequency of 2 ps⁻¹. The particle mesh Ewald method (Berendsen et al., 1984) was used for handling long-range interactions with a cutoff of 8.0 Å. The SHAKE algorithm (Ryckaert et al., 1977) was used to constrain bonds to hydrogen. The simulations were performed with an integration time step of 2 fs. Initially, we minimized the structures three times through 3,000 steps of steepest descent followed by 2,000 steps using the conjugate gradient algorithm. First, C_α atoms, the metal center, and the heavy atoms of the substrates were held fixed with a force constant of 5 kcal mol⁻¹ Å⁻². Second, only water molecules and C_α atoms were held fixed and, third, no restraints were used. Keeping the positions of the C_α atoms, the metal center, and the heavy atoms of the substrates fixed with a force constant of 2 kcal mol⁻¹ Å⁻², the systems were heated to 300 K over 25 ps of NVT-MD (constant number of particles, volume, and temperature) and simulated for further 25 ps. For density adaptation, 450  ps of NPT-MD (constant number of particles, pressure, and temperature) at 1.0  bar were performed. Afterward, the restraints on the atomic positions were gradually removed over 1,500 ns of NPT-MD. Finally, the systems were simulated for a further 1,000 ps without any restraints.

For production, the thermalized systems were simulated for a further 500 ns of NPT-MD without any restraints. For each system, five separate thermalizations and production simulations were performed, yielding an aggregate production simulation time per system of 2.5 ms. CPPTRAJ (Roe and Cheatham, 2013) was used to analyze the trajectories. Visualization of the molecular structures was done with PyMOL (DeLano, 2002).

CNA

Rigidity analysis was performed using the CNA software package (Pfleger et al., 2013). CNA is a front- and back-end to FIRST (Jacobs et al., 2001; Pfleger et al., 2013), which is used to construct networks of covalent and noncovalent interactions from input structures (Hermans et al., 2017). The strength of hydrogen bonds and ionic interactions is calculated by FIRST. Hydrophobic interactions are considered if the van der Waals spheres of two atoms are <0.25 Å apart. By successively removing the constraints below a certain energy cutoff E_cut (−7.0 kcal mol⁻¹ ≤ E_cut < −0.1 kcal mol⁻¹, ΔE_cut = 0.1 kcal mol⁻¹) from the network topologies, “constraint dilution trajectories” are generated (Radestock and Gohlke, 2008, 2011). The changes of the biomolecular stability along the “constraint dilution trajectories” are quantified in neighbor stability maps on a residue-wise level, yielding potential energy for each rc formed between pairs of residues (Rathi et al., 2015; Wifling et al., 2019). A residue’s chemical potential energy due to non-covalent bonding is obtained as the sum of all short-range rc the residue is involved in (Pfleger et al., 2013, 2017) (Eq. 1).

$$E_{i,\mathrm{CNA}}\mathrm{ }=\mathrm{ }\frac{1}{2}{\sum }_{j\ne i}^n{rc}_{ij,\mathrm{neighbor}}$$

(1)

The CNA software was applied on conformational ensembles of 2,500 snapshots extracted from the MD simulations of each system. Prior to the rigidity analysis, water molecules, sodium ions, and cap groups were removed from the structures, while NAD⁺, L-malate, and the Mg²⁺ ion were kept. Visualization of the molecular structures and the representation of the rc were done with PyMOL (DeLano, 2002).

Assessing the impact of individual substitutions in the GgNAD-MEβ1 isoform on the biomolecular rigidity and flexibility

To assess the impact of individual substitutions in the GgNAD-MEβ1 isoform on the biomolecular rigidity and flexibility, we applied a perturbation approach implemented in CNA (Pfleger et al., 2017). Ground state network topologies were generated from conformational ensembles as described above. Then, all covalent and noncovalent restraints associated with the side chain of a selected residue are removed, yielding the perturbed state, where the selected residue is stripped down to an alanine residue. As both states are generated from the same conformational ensemble, any changes in the local stability can be related directly to the in silico alanine substitution (Pfleger et al., 2017). Decomposing the free energy associated with the change in biomolecular stability on a per-residue level (ΔG_i,CNA) allows monitoring how the perturbation effect percolates through the enzyme structure, indicating the structural coupling of remote areas to the substitution site. By performing the same perturbation on corresponding residues in GgNAD-MEβ1 and GgNAD-MEβ2, ΔΔG_i,CNA is obtained for each residue i, which shows the changes in rigidity for the substitution between the isoforms (Eq. 2).

$${\Delta \Delta G}_{i,\mathrm{CNA}}\mathrm{ }=\mathrm{ }{\Delta G}_{i^{\mathrm{\text{'}}},\mathrm{ }\mathrm{CNA}}^{Gg\mathrm{NAD}-ME\beta 2}-\mathrm{ }{\Delta G}_{i,\mathrm{ }\mathrm{CNA}}^{Gg\mathrm{NAD}-ME\beta 1}$$

(2)

Here, i and i' are corresponding residues in GgNAD-MEβ1 and GgNAD-MEβ2 according to a sequence alignment (Tronconi et al., 2020).

Immunolocalization studies

Leaf tissue sampled from the middle of the youngest fully expanded leaf was fixed in 0.5% (w/v) glutaraldehyde, 4% paraformaldehyde (w/v) in sodium cacodylate buffer and prepared for transmission electron microscopy (TEM) and immunohistochemistry as described (Khoshravesh et al., 2017). Specific antibodies were raised against GgNAD-ME subunits based on small unique peptide sequences as antigens in rabbits. GgNAD-MEα and GgNAD-MEβ2 antibodies were produced by Agrisera (Vännäs, Sweden) in the project 2015-193. Anti-GgNAD-MEα antibody was produced with the antigen peptide “CPEYPTLVYKNA”, serum of ID 15D342 was used here. Anti-GgNAD-MEβ2 antibody was produced with the antigen peptide “DLTEFPGLVC”, serum ID 15D3304 was used. GgNAD-MEβ1 antibodies were produced by Cambridge Research Biochemicals (Billingham, UK) in the project CPVA17580 using “DLTEFPGLV” peptide as antigen, serum ID J30 was used here. The dilutions used were 1:100 (GgNAD-MEα), 1:50 (GgNAD-MEβ1, GgNAD-MEβ2), and 1:20 (secondary antibody; 18-nm colloidal gold-affiPure goat anti-rabbit IgG; Jackson Immunoresearch). Incubation times in the primary and secondary antibodies were 2 and 1 h, respectively. Controls were run by omitting primary antibody. TEM images were captured on a Phillips 201 transmission electron microscope equipped with an Advantage HR camera system (Advanced Microscopy Techniques, Woburn, MA, USA).

Statistical analyses

Statistical analysis for all experiments was performed using GraphPad Prism version 6.01 software using a two-way analysis of variance (with α < 0.05) and Tukey’s correction for post-hoc multiple comparisons. Results of multiple comparisons with confidence intervals and significance levels are shown in Supplemental Table S1.

Data availability

The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD026547 and 10.6019/PXD026547.

Accession numbers

Sequence data from this article can be found in the Arabidopsis Genome Initiative or GenBank/EMBL libraries under the following accession numbers: AtNAD-ME1 (α-subunit; At4g13560) and AtNAD-ME2 (β-subunit; At4g00570). The sequences of GgNAD-MEα (C. gynandra_Ggy11557), GgNAD-MEβ1 (C. gynandra_Ggy18870), GgNAD-MEβ2 (C. gynandra_Ggy19628), CaNAD-MEα (C. angustifolia_TRIN_DN31760c0g1i), CaNAD-MEβ1 (C. angustifolia_TRIN_DN18166-23437-30878-15844c0g1i), CaNAD-MEβ2 (C. angustifolia_TRIN_DN31681c1g5i), ThNAD-MEα (T. hassleriana_Th12536), ThNAD-MEβ1 (T. hassleriana_Th03046), and ThNAD-MEβ2 (T. hassleriana_Th09126) were acquired from transcriptome data, verified and correctly assembled using cDNA-based sequencing as the transcriptomes showed misassembled transcripts for several NAD-ME genes (Tronconi et al., 2020).

Supplemental data

The following materials are available in the online version of this article.

Supplemental Figure S1. Behavior of NAD-ME entities of Cleome species.

Supplemental Figure S2. Representative SDS–PAGEs showing purification steps of recombinant NAD-ME from E. coli BL21 cells transformed with different plasmids.

Supplemental Figure S3. Coomassie-stained native PAGE of purified recombinant GgNAD-MEs.

Supplemental Figure S4. Representative SDS–PAGEs showing different purification steps of the recombinant NAD-ME single proteins of T. hassleriana and G. gynandra.

Supplemental Figure S5. Structural stability of the simulated systems.

Supplemental Figure S6. Constraint network analysis.

Supplemental Figure S7. Impact of the substitutions in GgNAD-MEβ₁ on the structural rigidity.

Supplemental Figure S8. Transmission electron micrographs of BSCs of G. gynandra showing mitochondria with immunogold labeling.

Supplemental Table S1. Summary of the analysis of variance and Tukey’s multiple comparisons test for data in Figure 1.

Supplemental Table S2. Identification of NAD(P)-ME subunits by MS.

Supplemental Table S3. Kinetic parameters of recombinant GgNAD-ME entities in the absence and the presence of aspartate.

Supplemental Table S4. List of primers used in this study.

Supplemental Data Set S1. Influence of selected metabolites on the activity of GgNAD-ME isoforms.

V.G.M. and M.A.T. designed and supervised the research. M.H., M.A.T., J.P.S., T.L.S., M.S., G.P., and D.B. performed the research. G.P. performed the data curation. H.G. and D.B. performed the formal analysis. All authors were involved in writing the manuscript. V.G.M. and T.L.S. acquired funding.

The author responsible for distribution of materials integral to the findings presented in this article in accordance with the policy described in the Instructions for Authors (https://academic.oup.com/plcell) is: Veronica G. Maurino ().

Acknowledgments

H.G. is grateful for computational support by the “Zentrum für Informations und Medientechnologie” at the Heinrich-Heine-Universität Düsseldorf and the computing time provided by the John von Neumann Institute for Computing (NIC) on the supercomputer JUWELS at Jülich Supercomputing Centre (JSC) (user IDs: HKF7, VSK33).

Funding

This work was supported by grants of the Deutsche Forschungsgemeinschaft (MA2379/18-1 and the Germany’s Excellence Strategy EXC1028) to V.G.M., and the Canadian Natural Science and Engineering Research Council (2015-04878) to T.L.S.

Conflict of interest statement. The authors have no competing interests (financial/nonfinancial) that might be perceived to influence the interpretation of the article.

Author notes

†

Senior author.

These authors contributed equally (M.H., M.A.T.).