Stress-Induced Detoxification Enzymes in Rice Have Broad Substrate Affinity.

Reactive carbonyl compounds (RCCs) such as hydroxynonenol, malondialdehyde, acrolein, crotonaldehyde, methylglyoxal, and glyoxal accumulate at higher levels under stress in plants and damage the cell metabolic activities. Plants have evolved several detoxifying enzymes such as aldo-keto reductases (AKRs), aldehyde/alcohol dehydrogenases (ALDH/ADH), and glyoxalases. We report the phylogenetic relationship of these proteins and in silico analysis of rice-detoxifying protein structures and their substrate affinity with cofactors using docking and molecular simulation studies. Molecular simulations with nicotinamide adenine dinucleotide phosphate or glutathione cofactor docking with commonly known reactive substrates suggests that the AKRs, ALDH, and ADH proteins attain maximum conformational changes, whereas glyoxalase has fewer conformational changes with cofactor binding. Several AKRs showed a significant binding affinity with many RCCs. The rice microarray studies showed enhanced expression of many AKRs in resistant genotypes, which also showed higher affinity to RCCs, signifying their importance in managing carbonyl stress. The higher expression of AKRs is regulated by stress-responsive transcription factors (TFs) as we identified stress-specific cis-elements in their promoters. The study reports the stress-responsive nature of AKRs, their regulatory TFs, and their best RCC targets, which may be used for crop improvement programs.

Introduction

When plants are exposed to diverse environmental stresses, many metabolic, toxic intermediate compounds accumulate due to oxidative stress and inefficiency of catabolic enzymes. Some of these intermediate compounds are highly reactive, and they diffuse across cell organelles with a longer half-life and modify proteins.[1,2] Oxidative stress is ubiquitous, and the reactive oxygen species (ROS) damage proteins, lipids, and DNA, which leads to the production of reactive carbonyl compounds (RCCs). Several RCCs such as malondialdehyde (MDA), 4-hydroxy nonenal (HNE), acrolein, 4-oxo trans 2-nonenone (ONE), crotonaldehyde, glyoxal, glucosone, and deoxy glucosone accumulate under stress conditions.[2,3] These compounds possess highly electrophilic C=O groups. Furthermore, these compounds react with proteins and form protein carbonyls (PCs) that inactivate the function of metabolic enzymes and damage cellular homeostasis.[4] The RCCs carbonylate the proteins and inactivate their function by forming an adduct. Recent studies have shown that several key proteins involved in photosynthesis, Calvin cycle, anthocyanin biosynthesis, antioxidant mechanisms, protein synthesis turnover, etc. are carbonylated under stress conditions.[2] The rice genotypes exposed to accelerated aging showed higher accumulation of RCCs and loss of seed viability and seedling vigor.[5] Therefore, the carbonyl stress caused by RCCs is the major contributing factor for the diminished cell metabolism under stress conditions, thus affecting crop growth and productivity. From this context, it is inferred these compounds need to be scavenged to achieve cellular tolerance and improve resilience to stress conditions.

Plants have evolved many scavenging strategies to detoxify these cytotoxic compounds using enzymatic and nonenzymatic methods.[6,7] The antioxidants such as glutathione (GSH), ascorbic acid, and tocopherol, have shown potential in scavenging these RCCs.[8,9] Recently, several small molecules have been identified, which have the potential to detoxify the RCCs.[10,11] Apart from these, the cofactor-dependent enzyme families like nicotinamide adenine dinucleotide phosphate (NADPH)-dependent AKRs, aldehyde dehydrogenases, and GSH-dependent glyoxalases have shown a potential impact on enzymatic detoxification of the RCCs.[9,12,13]

The AKRs bind to RCCs like aldehydes and ketones and convert them into simpler alcohols using NADPH as the cofactor. Similarly, alcohol dehydrogenases bind to specific aldehydes and convert them into alcohol. Alcohol dehydrogenases convert higher alcohols into aldehydes and subsequently to less toxic alcohols. The glyoxalases convert MG to lactaldehyde and glyoxal using GSH as a cofactor.

The relevance of RCC-detoxifying enzymes has been well documented. The overexpression of AKR1 in rice has shown detoxification of RCCs, glyphosate, and NaCl-induced carbonyl compounds in rice and tobacco.[5,6] Overexpression of AKR1 showed detoxification of MG in rice suspension cells and tobacco and also showed tolerance to methyl viologen (MV)-induced oxidative stress.[14] Overexpression of ALDH7 in rice showed detoxification of RCCs and improved the seed viability.[15] Overexpression of glyoxalase I and glyoxalase II showed improved tolerance in tobacco for NaCl-induced RCCs by reducing MG.[13,16]

The diverse glycoxidation- and lipoxidation-induced RCCs have different reactivities to the target molecules. Though AKRs and other detoxifying enzymes have broad-spectrum substrate specificity, these enzymes may differ in their affinity to these diverse RCCs. Rice has 27 AKR family genes and glyoxalases and also other RCC-detoxifying enzymes. Information about their expression pattern, substrate specificity, and efficiency has not been studied. The broad-spectrum substrate specificity has made it difficult to identify candidate AKRs that can be targeted for breeding or genetic engineering. Because several diverse RCCs are generated under stress, a coordinated expression of several of the AKRs and other detoxifying enzymes may be necessary. From this context, identifying the relevant transcription factor (TF) that regulates many AKRs has significance because they can regulate a large number of downstream target genes.[17−19] Overexpression of TFs leads to upregulation of several antioxidative enzymes and RCC-detoxification genes, resulting in improved resistance.[6,19]

The oxidative-mediated changes in the gene expression of many regulatory genes occur at the very early stages of stress and they bind to specific clusters such as reactive oxygen species element (ROSE) cis-elements. The regulons related to ROSE binding sites have been identified in many promoters of genes that are upregulated under oxidative stress.[20] Several families of TFs, such as MYB, WRKY, zinc transporter, heat shock transcription factor, and basic-region leucine zipper (bZIP), have been found to be involved in the regulation of transcription of many genes, in response to oxidative stress.[21,22] Similarly enhanced oxidative stress tolerance is reported in bZIP activation in Brachypodium distachyon,[23] and overexpression of StDREB1[10] and OsGRAS23 results in salt and drought stress tolerance.[24]

The carbonyl stress is induced by a wide range of RCCs generated during the glycoxidation process. In addition, it is evident that plants also evolved several RCC-detoxifying mechanisms mainly by overexpressing several AKRs and other detoxifying enzymes. However, the substrate specificity of these enzymes and their affinity to RCC is not clearly elucidated so far. In spite of the broad substrate specificity of these enzymes, coordinated expression of these enzymes is crucial to achieving stress-induced oxidative stress damage in plants. For this, we provide information on the in silico structural variations among these groups of proteins. Our study indicates that they have a different level of affinity to different RCCs and have the potential to scavenge. In addition, our study provides leads in identifying the TFs that regulate AKRs based on the promoter analysis for the presence of stress-specific cis-elements. The expression of TFs correlates to the expression of AKRs in contrasting rice genotypes.

Results

RCC-Detoxifying Proteins Have a Similar Sequence Homology

Rice has 27 aldo–keto reductases (AKRs)[14] along with alcohol and aldehyde dehydrogenases (ADH/ALDH). In addition to these enzymes, the glyoxalase family of detoxifying enzymes is also involved in the detoxification of MG, an RCC. To know the sequence homology among AKRs and differences with ADH/ALDH and glyoxalases, analysis of the amino acid sequence homology was carried out using MEGA6. The sequence alignment of AKRs indicates highly conserved regions. The phylogenetic analysis clearly distinguishes the AKRs from other group of enzymes. Glyoxalases form a different group in the nearest neighbor analysis with a bootstrap value of 0.055. The aldose reductase has a close homology with the AKR Os05g0496200; otherwise, ADH/ALDH also forms a different group with a bootstrap value of 0.051 from AKRs (Figure A). The phylogenetic clade analysis also suggests that AKRs form a unique clade compared to ADH and ALDH family members. However, it is interesting to observe that a cinnamyl alcohol dehydrogenase forms a separate clade with glyoxalases (Figure B). These results demonstrated that all three groups of detoxifying enzymes have unique structural features and AKRs are unique and may bind effectively to RCCs compared to other proteins.

Figure 1

Phylogenetic analysis of rice RCC-detoxifying enzymes. (a) Phylogenetic analysis of AKRs, ALDH, ADH, and glyoxalases using amino acid sequences from the rice database by the nearest neighbor joining method using MEGA6 with 1000 bootstrap values and (b) phylogenetic clade analysis indicating the different families of proteins differing based on their amino acid homology.

Detoxifying Enzymes Attain Conformational Changes upon Cofactor Interaction

To know the structural variations in RCC-detoxification enzymes, the 3D protein structures were predicted for AKRs, ALDH, ADH, and glyoxalase using the I-TASSER web tool. As predicted by the protein database, all the AKRs possessed an eight α–β barrel architecture (Figure S1A). Similarly, ADH and ALDH also possessed the α–β barrel motif (Figure S1B,C). However, the glyoxalases possessed the βαβββ architecture, which is a characteristic feature of these groups of proteins (Figure S1D).

The cofactor binding to these enzymes resulted in conformational changes. The NADPH binding pocket in AKRs is localized within the carboxy-terminal phase of the central β-barrel (Figure A). The cofactor binding affinity was estimated using protein–NADPH docking. Different amino acids in each AKR are found to interact with NADPH to form an apoprotein and a holoenzyme. Among the AKRs, Os03g0237100 has the highest affinity with a docking score of −11.022 and interacts at the amino acid residues R194, R197, R201, and V205. The AKR Os02g0123500 has the lowest docking score of −6.9442 and interacts at the amino acids E116, N145, K161, H180, T181, and G182 (Table S1A). These observations suggest that different amino acids in the AKRs participate in cofactor NADPH binding with different affinity levels (Table and Table S1A).

Figure 2

Cofactor binding in RCC-detoxifying enzymes. (A,B) Cofactor NADPH binding in the center of β barrel in AKRs and ADH proteins, (C) NADPH binding at the proximal end of the protein in ALDH, and (D) cofactor GSH binding site in glyoxalase I at the proximal end of the protein. The protein structures were built using RaptorX. The docking studies were performed using Glide Schrödinger 2017–3 (extra precision docking).

Table 1

Affinity of AKRs with their Respective Cofactor and RCCs

	detoxifying gene	affinity to cofactor/binding score	affinity to RCC	binding pockets
1	Os05g0474600	Medium/ −8.926	Medium	Medium
2	Os02g0817500	High/ −10.462	High	High
3	Os04g0338000	Medium/ −9.666	Medium	High
4	Os04g0447700	Medium/ −9.029	High	Low
5	Os10g0113000	Low/ −7.7693	Low	Medium
6	Os10g0113900	High/ −10.265	Medium	High
7	Os10g0419100	High/ −10.251	Medium	Medium
8	AK073738	Low/ −7.581	Medium	Medium
9	Os03g0237100	High/ −11.022	Low	High
10	aldehyde dehydrogenase 7 Os09g0440300	Low/ −5.787	High	Medium
11	alcohol dehydrogenase 1 Os11g0210300	Medium/ −9.084	Medium	High
12	glyoxalase I Os08g0191700	Low/ −7.614	Medium	Medium

ADH has a central NAD binding pocket similar to AKRs (Figure B). However, ALDH showed a cofactor NADP binding pocket at the proximal side (Figure C). Among the ADH/ALDH tested, cinnamyl alcohol dehydrogenases showed interaction with NADP at N166 and T302 amino acids with a docking score of −10.385. Similarly, ALDH7 has NADP binding sites at R129 with a docking score of −5.785, which is less in these groups of proteins (Table S1B). The glyoxalases showed a cofactor GSH binding pocket at the proximal end (Figure D). The Gly I showed GSH interactions at amino acids R27, N28, F79, E157, K207, and N208 with a docking score of −7.614 and Gly II showed cofactor interaction at I49, V54, R71, L73, and N75 amino acids with a docking score of −7.440 (Table S1C).

The cofactor binding to apoprotein brings in conformational changes, which makes the active sites accessible for the substrate to bind and interact for catalysis. AKRs have many active sites where diverse substrates can bind, for example, AKR1 has nine active sites that expand upon cofactor interaction, and they were identified at different structural depths (Figure S2 (i)). The expansion of enzyme active sites in terms of depths was higher in the case of ALDH7 compared to AKR, ADH, and glyoxalases (Figure S2 (ii,iii)). There were minimal structural changes in the active sites in glyoxalase I with cofactor GSH (Figure S2 (iv)).

The cofactor interaction and conformational changes were reconfirmed by root-mean-square deviation (RMSD) using molecular dynamics (MD) simulations. The RMSD analysis shows that ALDH7 attains the maximum conformational changes compared to AKRs and glyoxalases (Figure ). However, the GSH binding to glyoxalase I did not yield significant structural changes (Figure ). The structural analysis and MD simulations indicate that substrate-binding pockets for AKRs are relatively wide with a deep elliptical cavity. Several AKRs exhibit many of these cavities, indicating that they could bind to different substrates at different locations.

Figure 3

RMSD plot for the Cα atom depicting the structural changes on binding with the cofactor for 10 ns with different proteins. The simulation studies were carried out using the Nanoscale Molecular Dynamics program (NAMD) (Linux-x86_64 multicore).

AKR–NADPH Complexes Interact with RCCs with Different Affinities

The RCCs such as MDA, HNE, ONE, glyoxal, MG, acrolein, crotonaldehyde, glycolaldehyde, glucosone, 3-deoxy glucosone, and 3-deoxy fructose have shown to accumulate at higher levels under stress conditions and react with proteins and other macromolecules forming PCs, aggregates, advanced lipoxidation end products, and advanced glycation end products. To detoxify these RCCs, the cofactor binding to the scavenging enzymes brings in conformational changes that allow substrates to bind at the active sites.

To assess the affinities of diverse AKRs, ADH, ALDH, and glyoxalases, the apoprotein cofactor (NADPH/NAD/NADP/GSH) was docked with 11 most commonly found RCCs using the Schrödinger tool. All the AKRs were found to interact with all the RCC substrates at different amino acids (Figure A and Table ). Each RCC substrate has a different binding location on AKR–NADPH complexes (Table S2A). In the case of AKR1, MDA binds to only one amino acid R286 with a docking score of −4.5, whereas HNE could interact with nine amino acids in AKR1 with a docking score of −4.7. AKR1 interacts with ONE at N156, Y203, G290, and W307 amino acids and the interactions were found to be stable when analyzed using MD simulations for 2 ns. AKR1 interacts with crotonaldehyde at Trp307 with a sigma bond at the first carbon position and at Arg 283 with a hydrogen bond (Figure S3 (i)). The docking score of AKR1, to different RCCs, is in the range of −3.7 to −5.2 with ONE and crotonaldehyde, respectively (Table S3A). Among the AKRs, Os03g0237100 has the highest docking score of −6.2 with ONE. The least is for AKR Os04g0337500 with a docking score of −4.2 with deoxy-glucose (Figure A and Table S3A).

Figure 4

Protein–substrate interaction. (A) Docking of the AKR–NADPH complex with the substrate; the phylogenetic and graphical interactions were developed using R studio 1.1.456 (https://rstudio.com/products/rstudio/download/) and ggplot libraries (https://rstudio.com/products/rpackages/). (B) RMSD plot for the Cα atom, depicting the structural change on binding with respective RCCs with different proteins. Interaction of AKR1 with crotonaldehyde, ADH1 with 4-hydroxy-trans-2-nonenal, ALDH with 4-hydroxy-trans-2-nonenal, and glyoxalase with 4-Hydroxy-trans-2-nonenal.

The ADH interacts with HNE at Phe 149 and 67 with an alkyl bond with the first carbon atom and Trp54 at the hydroxyl group (Figure S3 (ii)) and at other amino acids H84, P131, T136, F157, and V299 with a docking score of −6 (Table S3B). ALDH7 also showed the highest docking score of −9 when HNE interacted at Phe 135 and cys 44 with an alkyl bond and formed a hydrogen bond with the OH group with Gly134 and formed a C–H bond at Gly68 (Figure S3 (iii)). The RMSD analysis shows that ALDH7 attains the maximum conformational changes compared to AKRs and glyoxalases (Figure B). The glyoxalases were found to interact with all the RCCs; however, their interaction with MG is less compared to HNE. GlyI showed hydrogen bonding with MG at Trp92 and has a docking score of −4.4 and Gly II showed a docking score of −3.5 (Figure S3 (iv) and Table S3C). The glyoxalase I showed a low RMSD with 4-HNE (Figure B).

MDA has the highest affinity toward ALDH7 with a docking score of −5.5 and the least is AKR Os04g0337500 with a docking score of −2.9. HNE, MG, and acrolein have the highest docking score among the AKRs with OsAKR1 with values of −4.7, −4.6, and −4, respectively. The AKR Os04g0337500 has docking scores of −4.2, −2.9, and −2.7 with HNE, MG, and acrolein, respectively. It is interesting to observe that all the AKRs have different levels of affinity to the RCCs tested (Figure A). Another interesting factor is that glyoxalases also showed a broad substrate affinity other than MG (Table S3C).

Figure 5

Interaction of AKR proteins with RCCs. (A) Network of molecules and proteins showing multiple substrate specificity and (B) specificity of protein–substrate interaction based on the highest docking score. Interactions of AKR proteins with RCCs were plotted using Cytoscape v3.6.1. The RCC compounds are present in the inner ring with AKRs on the peripheral ring.

Based on the docking score, many RCCs showed the highest affinity toward ALDH7 compared to both AKRs and glyoxalases (Table ). As the cofactor binding affinity exposes multiple binding pockets, the apoprotein may have the highest affinity toward substrates. ADH1 and AKRs Os02g0817500 and Os04g0338000 have medium affinity to RCC but they have more number of substrate-binding pockets. The AKRs Os10g0113900 and Os03g0237100 showed low affinity to RCC but have high binding pockets (Table ). These enzymes may work efficiently in an in vitro system; however, more studies are required to evaluate this hypothesis. Among the substrates, 4-HNE, ONE, glucosone, and 3-deoxy-glucosamine have the highest binding affinity to AKRs. Many proteins including AKRs, ALDH, and ADH have the highest affinity toward HNE and some proteins such as OsAKR1 have affinity toward crotonaldehyde (Figure B). 3-Deoxyfructose, MDA, and crotonaldehyde have medium affinity, and glyoxal, MG, acrolein, and glycolaldehyde have the lowest affinity toward many AKRs (Table S4). The in vivo effect of all these enzymes depends on their level of accumulation, cofactor availability, and site of their presence with the substrate.

Differential Expression of AKRs in Rice Is Regulated by Transcription Factors

To study the differential stress response of AKRs, 10-day-old seedlings from contrasting rice genotypes AC39020 (resistant) and BPT5204 (sensitive) were exposed to oxidative stress and the transcriptome profile was developed using the Agilent Microarray Platform (Genotyping technologies, https://www.genotypic.co.in). Among the several genes that were differentially expressed between the contrasting genotypes, the expressions of all AKRs and few ADH/ALDH and glyoxalases were assessed (Figure A). The transcript levels of 15 AKRs were higher in resistant genotype AC39020 and only seven AKR transcripts were more in the susceptible genotype BPT5204 (Figure A). The transcript levels of Os05g074600, Os07g0142900, and Os05g0456300 were found to be more than threefold, and other AKR genes showed more than onefold accumulation in the AC39020 resistant rice genotype. In the sensitive genotype, only Os07g014300 and Os07g0142900 showed more than twofold higher transcripts and the remaining genes showed less than onefold accumulation. A few AKR transcripts were reduced in both the genotypes, and among them, more number of transcripts were lower in the susceptible genotype BPT5204 compared to the resistant AC39020 genotype. Furthermore, the ALDH and ADH gene transcripts were also reduced in the sensitive genotype. However, ALDH, Os04g0447700, and Os02g0817500 transcripts showed about onefold accumulation in the resistant AC39020 genotype, whereas 0.5-fold higher accumulation was observed in BPT5204. The transcripts of glyoxalases were also higher in both the genotypes; however, the levels were more in the resistant AC39020 genotype (Figure A).

Figure 6

Expression of AKRs and TFs in oxidative stress. (A) Differential expression of AKRs in two contrasting rice genotypes differing in oxidative stress, that is, resistant AC39020 and sensitive BPT5204. (B) TFs binding to promoters of AKRs, ADH, ALDH, and glyoxalases. Promoters of AKRs with TF interaction—the 2 kb upstream sequences of all the RCC-detoxifying enzyme encoding genes were assessed using PlantPAN2.0. Based on the respective cis-elements present in the promoters, the arrows were drawn manually. (C) Differential expression of TFs in contrasting rice genotypes. The microarray study was carried out for 10-day-old seedlings exposed to oxidative stress, and differential expressions of TFs and AKRs were filtered for comparison over the control conditions.

Since the expression of AKRs varied in both rice genotypes under stress, we speculate that the cis-elements in these gene promoters are responsible for the differential expression. The 2 kb region of the upstream promoter sequence analysis of AKRs suggests that all AKRs have a few stress-responsive cis-elements, and especially, they are rich in bZIP, WRKY, and ERF family cis-elements. The TFs bZIP37, bZIP12, bZIP23, WRKY23, ERF53, and ERF54 have shown to bind to these cis-elements (Figure B). It is interesting to note that the OsbZIP23 binding cis-elements are present in almost all the AKRs (Table and Table S5). Furthermore, the expression of these six TFs in contrasting genotypes was assessed.

Table 2

Expression of Detoxifying Genes in Contrasting Rice Genotypes under Stress and Presence of Specific Cis-Elements in Their Promoters

	detoxifying gene	TF binding/upregulation in rice genotype	AC39020	p value	BPT5204	p value
1	Os05g0474600	bZIP12 and 23/AC39020	2.46	0.00013	0.83	0.002175
2	aldehyde dehydrogenase 7 Os09g0440300	WRKY23 and bZIP23/AC39020	1.58	0.01370	–0.4	0.0007937
3	alcohol dehydrogenase 1 Os11g0210300	WRKY 23 and bZIP12 and 23/AC39020	1.41	0.00004	0.46	0.0073863
4	glyoxalase I Os08g0191700	WRKY 23 and bZIP12 and 23/AC39020	1.16	0.004556	0.89	0.0008192
5	Os02g0817500	WRKY 23 and bZIP12 and 23/AC39020	0.87	0.013857	0.39	0.0211906
6	Os04g0338000	WRKY 23 and bZIP12 and 23/AC39020	0.64	0.00014	–0.12	0.0353247
7	Os04g0447700	bZIP37 and ERF53 and 54/BPT5204	0.32	0.050511	0.47	0.0144040
8	Os10g0113000	bZIP37 and ERF53 and 54/BPT5204	–0.41	0.00259	–0.32	0.0280583
9	Os10g0113900	WRKY 23 and bZIP12 and 23/AC39020	–0.71	0.00562	–1.06	0.0014831
10	Os10g0419100	bZIP37 and ERF53 and 54/BPT5204	–1.04	0.018267	–0.68	0.0102426
11	AK073738	bZIP37 and ERF53 and 54/BPT5204	–1.98	0.009285	–0.62	0.0079127
12	Os03g0237100	bZIP37/BPT5204	–1.98	0.009285	–0.62	0.0079127

In the oxidative stress-induced microarray data from contrasting rice genotypes, the resistant genotype AC39020 showed higher transcript levels of TFs OsbZIP12 (Os01t086730), OsbZIP23 (Os02t076670), and OsWRKY23 (Os01t063700). The sensitive rice genotype BPT5204 showed higher levels of transcripts of TFs OsbZIP37 (Os04t063700) and ERF53 and 54 (Os01t065740 and Os01t022410) (Figure C). The data suggest that in resistant genotype AC39020, the number of AKRs that were upregulated at a higher fold could be due to the upregulation of these TFs and some of the AKRs that were upregulated in the sensitive genotype may be correlated to the upregulation of TFs. However, this upregulation was not sufficient to enhance the transcript levels of any AKRs in the sensitive genotype as observed from the microarray data (Table ).

Discussion

During stress conditions, the levels of ROS and reactive aldehydes accumulate and cause damage to proteins. Thus, introduction of antioxidant mechanisms or detoxifying systems is essential for higher levels of tolerance by scavenging ROS or detoxifying the reactive aldehydes. The reactive compounds such as HNE, acrolein, MDA, and MG are known to accumulate both in mammals as well as in plants, which cause modifications in nucleic acids and proteins and leads to the inactivation of proteins.[2] Plants have evolved several adaptive mechanisms to survive and develop defense strategies. Both biotic and abiotic factors disrupt the delicate balance between endogenous levels of ROS and oxidants.[25,26] There are three different enzyme systems such as long-chain ALDH/ADH, short-chain dehydrogenases/reductases, and AKRs. Interestingly, all the proteins possess Rossmann fold for NADPH cofactor binding to reduce RCCs. The presence of α–β barrel motifs in the 3D structure is a characteristic feature of all these three groups of proteins. However, glyoxalases have α,-β,β,β motifs. Among them, AKRs are a huge family in rice, and they are known to reduce aldehydes and ketones. Apart from these, they also rarely reduce monosaccharides, steroids, prostaglandins, and polycyclic hydrocarbons.[27,28] The bioinformatic analysis of AKRs, ALDH, and glyoxalases suggests that all AKRs form a unique clade (Figure ). Similar observations have been reported from Arabidopsis and other species.[28] All these enzymes have a broad physiological activity and substrate specificity; however, it is not clear why plants have so many AKRs and other detoxifying enzymes.

The substrate specificity of AKRs including plants and mammals suggested that these groups of enzymes bind to the compounds which contain C=O reactive aldehydes and ketones. The lipid peroxidation compounds like MDH, HNE, and acrolein and glycation products like MG, glucosone, etc. have these carbonyl atoms and they possess high reactivity for proteins.[2,28] The AKRs are NADPH-dependent and use cofactors for catalysis. Based on the active site, or the loops created by NADPH binding, the conformational changes determine the substrate-binding efficiency. As we observed by protein–NADPH docking, several amino acids from each AKR possessed NADPH binding ability (Table ). Similarly, cofactor binding for different AKRs at different amino acid locations was reported.[29] Upon cofactor binding, ALDH7 attained maximum conformational changes compared to glyoxalases with GSH (Figure ). Structural changes occur as a result of cofactor binding, which is evidenced through many X-ray structure measurements from members of AKRs.[29] Similar structural changes were observed in bioinformatic docking studies with apoenzymes (Figure ). However, the overall protein backbone did not change very much with the root-mean-square (RMS) difference between the Cα atoms. The RMS values for glyoxalases are significantly lower, indicating that they have specific substrates, whereas other detoxifying proteins with higher RMS values can bind to multiple substrates. It has been observed that upon cofactor docking, some residues of side chains undergo minor movements. Similar observations were noticed in the case of AKR1B.[30] The catalytic mechanism of AKRs follows an ordered bi-bi reaction, which is the hallmark feature of all AKRs and is conserved.[31]

The in silico docking of the NADPH apoenzyme with RCC compounds suggests a different set of amino acids interacting with substrates. However, the binding efficiency or affinity to substrates varies with differential docking scores. Based on our analysis, we could identify the best substrate for each AKR and other detoxifying proteins, which are highly essential to know the priorities of AKRs. ALDH7 showed the highest affinity to many RCC substrates and HNE was found to have the highest binding score for many enzymes (Table ). However, the data indicate that all the AKRs possessed affinity for all the substrates tested but the level of affinity varied with substrates. The RCCs are known to accumulate under stress conditions; as a result, several AKR genes in rice have been upregulated along with other differential responsive genes.[14] However, it is interesting to mention that MG, a known substrate for glyoxalase with GSH, also showed some affinity toward other carbonyl-containing molecules. Further in vitro studies are needed to confirm the substrate specificity of this glyoxalase against other RCCs.

The oxidative stress-tolerant rice genotype AC39020 may have a better detoxification ability compared to the sensitive BPT5204 genotype. The resistant genotype should have more number of genes expressed and their expression could be under stress-induced TFs. The AKRs were found to be stress-responsive as evident from the microarray data (Figure A). The resistant genotype AC39020 and 17 detoxifying gene transcripts were found to accumulate more than twofold, and in sensitive BPT5204 genotype, only 12 genes had higher transcripts. The regulation of the AKR gene expression may be controlled by upstream TFs.[18] The promoter analysis of AKR genes has identified the stress-specific cis-elements. The cis-elements for TFs such as OsbZIP37, OsbZIP12, WRKY23, OsbZIP23, ERF 53, and ERF 54 are found to be predominately present in the promoters of many detoxifying genes. The TFs bZIP12, bZIP23, and WRKY23 were upregulated in resistant genotype AC39020. The overexpression of OsbZIP12 in rice shows improved tolerance to drought, and seedlings are hypersensitive to ABA.[32] OsbZIP23 transgenics are sensitive to ABA and they confer resistance to salinity and drought in rice.[33] Similar heterologous expression of OsWRKY23 in Arabidopsis showed dark-induced senescence and pathogen defense.[34] The overexpression of TFs upregulates many downstream target genes.[17] Our data for resistance genotype correlate with the expression of TFs and AKRs suggesting that AKRs could play a major function in improving the resistance to oxidative stress tolerance. The TFs bZIP37, ERF53, and ERF54 have shown cis-elements that are present in the promoters of few AKRs that were upregulated in the sensitive genotype BPT5204; however, the number of AKRs upregulated is very less. The expression of bZIP37 (OsTGAP1) causes elicitor-induced hyperaccumulation of diterpenoid phytoalexins in rice cells, which showed a defense response.[35] The role of ethylene-responsive TFs (ERF) in regulating redox genes, ABA, and Jasmonate signaling genes has been well reported.[36] In the resistant AC39020 rice genotype, the higher expression of many AKRs, ALDH7, ADH1, and glyoxalase I is due to the presence of cis-elements of bZIP12, bZIP23, and WRKY23 (Table ). These TFs also showed significant upregulation under oxidative stress. The expression of TFs is tightly correlated with the expression of AKRs in the respective genotypes. However, in the resistant genotype, the TFs binding to the promoter may have resulted in higher expression of AKR transcripts. These highly expressing AKRs showed higher conformational changes with cofactor NADPH and also showed higher affinity to many RCCs (Tables and 2), indicating that they could efficiently scavenge these cytotoxic compounds and provide improved resistance in genotype AC39020.

Our study demonstrates that AKRs have broad substrate specificity and they are stress-responsive. The stress-responsive TFs induce the AKRs that have high conformational ability and high affinity toward RCCs and hence the resistance. The TFs play a crucial role in regulating the detoxification genes and RCCs. These genes could be potential molecular markers for crop improvement programs.

Materials and Methods

Protein Structure Prediction, Molecular Docking, and Simulations

The protein sequences for the AKRs were obtained from the National Centre for Biotechnology Information (NCBI) protein sequence database and substantiated at UniProt. The protein sequence was utilized to build the secondary structure using I-Tasser (https://zhanglab.ccmb.med.umich.edu/I-TASSER/). The final model chosen was evaluated before and after refinement. The stereochemical properties were calculated and tested for the structure quality using the PROCHECK and Verify3D environment profile analysis methods before submitting to the Structural Analysis and Verification Server (SAVES v5.0) (https://servicesn.mbi.ucla.edu/SAVES/). Based on the results, the protein structure was refined to obtain reliable scores with residues in the allowed region >90% and 2D-to-3D conversion rate > = 0.8 on a scale of 0 to 1. The Schrödinger platform, Small Molecule Discovery Suite 2017–3, was used for the docking experiment. The glyphosate and PsAKR1 with NADPH were used as the model ligand–enzyme pair .[37] Various tools like Protein Preparation Wizard, LigPrep, SiteMap, Receptor grid generation, and Glide were used.[38−41]

The protein was subjected to preprocessing for defining residues and heteroatoms and for verifying valency. Under structure refinement, hydrogen bond assignment was carried out at neutral pH. Structure minimization of less than 0.30 Å was conducted with OPLS_2005 for a defined force field. The top ranked potential binding sites were selected with restricted hydrophobicity. LigPrep was used to generate 3D structures of ligand molecules with 32 conformers for each. An extra precision ligand docking was performed using Glide with a scaling factor of 0.8 and a partial cutoff of 0.15. Epik[42] and Impact state penalties were added to the docking score.

MD Simulations

MD simulations for the modeled proteins were performed using the program NAMD v 2.9 ,[43] and all files were generated using visual molecular dynamics.[44] The protein was solvated with a transferable intermolecular potential (TIP3P) solvent water box with a 5 Å layer of water for each direction of the coordinate structure. Molfacture was used to draw the structure of NADPH and GSH, and a force-field toolkit (ffTK)[45] was used to build the parameter file, geometry file, charge constraints, optimized bonding, and scan torsions. This step was crucial in containing the NADPH was part of simulation system. For the minimization and equilibration of proteins along with NADPH in the water box, we assumed force-field parameters excluding a scaling of 1.0 Å and a cutoff of Coulomb forces with a switching function.

MD simulations were performed by starting at 12 Å, reaching zero at a distance of 10 Å, and ending at 14 Å with a margin of 2.0 Å. Integrator parameters also included 2 fs/step for all rigid bonds and nonbonded frequencies were selected for 1 Å and full electrostatic evaluations for 2 Å were used, with 10 steps for each cycle. The particle mesh Ewald method was used for electrostatic interactions of protein system periodic boundary conditions with grid dimensions of 1.0 Å. To eliminate bad water constraints, the protein preliminary energy was minimized via 500 steps of the Powell algorithm. The temperature was reassigned after two frequency steps and was held constant at 310 K. Further simulations were carried out for 5,000,000 runs (a total of 10 ns) of Langevin dynamics to control the kinetic energy, temperature, and/or pressure of the system.

Network Analysis

The docking results of AKRs with RCCs were plotted as a heatmap using R software (https://www.R-project.org) with ggplot libraries.[46] Network analysis and differential gene expression of TFs were carried out to identify their binding to the promoters of AKRs, ADH, ALDH, and glyoxalases. The 2 kb upstream sequences of all the RCC-detoxifying genes were assessed using PlantPAN2.0 (http://plantpan2.itps.ncku.edu.tw/index.html). The sequences were added as an input to RSAT (http://rsat.eead.csic.es/plants/index.php) and Motif discovery–RSAT oligo analysis.[47] The obtained results were mapped back to the PlantPAN database to verify the presence of cis-elements. A network was obtained using Cytoscape.[48]

Promoter Analysis To Identify the TF Binding Sites in AKR Genes

The promoter sequences (2 kb upstream sequences) of all the AKRs were downloaded from the RiceXpro website for identification of stress-responsive cis-elements. These upstream sequences were analyzed for the presence of cis-elements using PlantPAN 2.0 (http://plantpan2.itps.ncku.edu.tw).

Imposition of Oxidative Stress to Contrasting Rice Genotypes

The contrasting rice genotype seedlings were grown in plastic bowls for 10 days in a Yoshida rice growth nutrient (YRGN) solution. Oxidative stress was imposed on the 10th day by replacing the medium with the YGRN solution supplemented with 0.2 μM MV + 50 mM NaCl + 50 μM CdNO3. Seedlings were allowed to grow under stress for another 8 days and sampling was done by collecting both the root and shoot together for mRNA isolation for microarray analysis. Total RNA from AC39020 and BPT5204 under control and stress conditions was isolated using a Spectrum Plant Total RNA kit (Sigma) and was quantified using a spectrophotometer (BioSpec-nano, Shimadzu). The integrity of total RNA was verified on an Agilent 2100 Bioanalyzer using the RNA 6000 Nano LabChip (Agilent Technologies).

Microarray Data Development and Analysis

The samples for gene expression were labeled using an Agilent Quick Amp labeling kit (p/n5190–0442). Each of the total RNA (500 ng) was reverse-transcribed at 40 °C using an oligodT primer with a T7 polymerase promoter and converted to double-stranded cDNA. The synthesized double-stranded cDNA was used as a template for cRNA generation. cRNA was generated by in vitro transcription and the dye Cy3 CTP (Agilent) was incorporated during this step. The cDNA synthesis and in vitro transcription steps were carried out at 40 °C. The labeled cRNA was cleaned up using Qiagen RNesay columns (Qiagen, Cat No: 74106) and quality-assessed for yields and specific activity using the NanoDrop ND-1000. Then, 500 ng of the labeled cRNA sample was fragmented at 60 ° C and hybridized on to an Agilent Custom Rice GXP 8×60K array designed by Genotypic Technology Private Limited.

The extracted raw data were analyzed using Agilent GeneSpring GX software. The rice whole genome 8×60K array covers chloroplast genes, mitochondrial genes, and coding region genes. Normalization of the data was carried out in GeneSpring GX using the 75th percentile shift. The genes were assigned to different pathways utilizing the KEGG database using the software “Genotypic Biointerpreter”—a biological analysis software. Biointerpreter is a user-friendly web-based biological interpretation tool developed by Genotypic Technology Private Limited, Bangalore.

The intensity values were log (log base 2)-transformed and compared between samples. A gene showing an increased expression of onefold (log base 2) and above under stress conditions compared to its value in control is considered as upregulated. Similarly, a gene showing a decreased expression by less than onefold is considered as downregulated. Student’s t-test was conducted among the replicates for statistical significance and a p value of 0.05 was fixed.

Differentially expressed genes (DEGs) between genotypes. The mean gene expression (averaged log2 values of two replications) values of AC39020 and BPT5204 were compared to get the DEGs. The mean log2 values of a particular gene were compared between two genotypes and the difference was calculated. The gene considered should have a p value less than or equal to 0.05. By default, if the log2 difference is more than 0.6, it is considered as differentially expressed. However, the stringency can be increased by considering a higher difference of 1.0 or more.