Methylglyoxal detoxifying gene families in tomato: Genome-wide identification, evolution, functional prediction, and transcript profiling

Methylglyoxal (MG) is a highly cytotoxic molecule produced in all biological systems, which could be converted into non-toxic D-lactate by an evolutionarily conserved glyoxalase pathway. Glutathione-dependent glyoxalase I (GLYI) and glyoxalase II (GLYII) are responsible for the detoxification of MG into D-lactate in sequential reactions, while DJ-1 domain containing glyoxalase III (GLYIII) catalyzes the same reaction in a single step without glutathione dependency. Afterwards, D-lactate dehydrogenase (D-LDH) converts D-lactate into pyruvate, a metabolically usable intermediate. In the study, a comprehensive genome-wide investigation has been performed in one of the important vegetable plants, tomato to identify 13 putative GLYI, 4 GLYII, 3 GLYIII (DJ-1), and 4 D-LDH genes. Expression pattern analysis using microarray data confirmed their ubiquitous presence in different tissues and developmental stages. Moreover, stress treatment of tomato seedlings and subsequent qRT-PCR demonstrated upregulation of SlGLYI-2, SlGLYI-3, SlGLYI-6A, SlGLYII-1A, SlGLYII-3B, SlDJ-1A, SlDLDH-1 and SlDLDH-4 in response to different abiotic stresses, whereas SlGLYI-6B, SlGLYII-1B, SlGLYII-3A, SlDJ-1D and SlDLDH-2 were downregulated. Expression data also revealed SlGLYII-1B, SlGLYI-1A, SlGLYI-2, SlDJ-1D, and SlDLDH-4 were upregulated in response to various pathogenic infections, indicating the role of MG detoxifying enzymes in both plant defence and stress modulation. The functional characterization of each of these members could lay the foundation for the development of stress and disease-resistant plants promoting sustainable agriculture and production.


Introduction
The effects of climate change due to global warming result in various forms of environmental stress, dramatically decreasing crop yield and declining economic development.During adverse conditions, cytotoxic methylglyoxal (MG) is produced in biological systems through several non-enzymatic reactions as a byproduct of carbohydrate, protein and fatty acid metabolism [1].MG introduces toxicity by generating free radicals and advanced glycation end-products (AGEs), thus reducing growth as well as productivity [2].Dicarbonyl methylglyoxal is a renowned mutagen and genotoxic reactive, capable of altering protein and nucleic acid.Any modification of arginine's guanidine groups in protein causes the production of AGEs, which cause MG-induced damage in biological systems [2].Imidazopurinone derivatives, the most common AGEs generated from MG, have been linked to genomic integrity loss and genotoxicity [3].In response to the herbicide imazethapyr, crops such as lentils accumulate methylglyoxal, detoxified by glyoxalase pathway enzymes [4].
The glyoxalase pathway, responsible for detoxifying MG, is highly conserved and ubiquitous among lower prokaryotes as well as higher eukaryotes such as humans and plants [5].Two conventional glyoxalase enzymes, glyoxalase I (GLYI) and glyoxalase II (GLYII), convert MG derivative hemithioacetal (HTA) into D-lactate through two sequential reactions in a glutathione-dependent pathway [6].By contrast, a unique glyoxalase enzyme, a homolog to DJ-1 proteins in Arabidopsis, glyoxalase III (GLYIII) catalyzes the whole process in a single step in a glutathione-independent manner [7][8][9] (Fig 1A).In the final step of MG detoxification, TCA cycle metabolite pyruvate is formed from D-lactate by the action of D-lactate dehydrogenase (D-LDH) [10].
The glyoxalase system is significant in maintaining the normal homeostasis of plants by preventing any alteration of molecular and morphological characteristics upon exposure to stress [6].It has been thoroughly investigated because of its capacity to protect against harmful environmental conditions by detoxifying MG, which would otherwise cause significant hazards under diverse abiotic and biotic stresses [18].The role of glyoxalase in transgenic plants confirmed its function in abiotic stress modulation.Overexpression of the GLYI or GLYII gene showed MG and salt tolerance in the transgenic tobacco [19,20].Transgenic tobacco lines overexpressing both GLYI and GLYII, on the other hand, showed greater tolerance to salt stress [19].Brassica juncea exhibited a considerable increase in GLYI expression against salt, mannitol, and heavy metal stress [19].Similarly, the GLYI members were substantially expressed in Manihot esculenta (cassava) in response to iron-induced stress, whereas the overexpression of MeGLYI-13 confers iron toxicity tolerance in engineered Arabidopsis [21].When exposed to several abiotic stressors, the rice GLYII gene showed early activation, including NaCl, evaporation, heat, and abscisic or salicylic acid treatments [22].Biotic stress related to Agrobacterium tumefaciens infectious altered glyoxalase expression patterns [23].GLYI genes had lower expression in rice after being infected with Xanthomonas oryzae pv.Oryzae or Pyricularia grisea [24].
Substantial research has been conducted on the glyoxalase pathway in various species of microbes.However, only a few plant species had been analyzed genome-wide to identify MG detoxifying gene family members.Tomato is an economically important crop cultivated worldwide.According to FAOSTAT of the United Nations, about 186,821,216 metric tonnes of tomatoes were produced worldwide in 2020 [25].However, the production suffered from exposure to adverse weather and pathogenic infection.Therefore, a genome-wide investigation has been conducted to evaluate the significance of glyoxalase and D-LDH members in the stress physiology of tomato, which had identified a total of 13 GLYI, 4 GLYII, 3 GLYIII, and 4 D-LDH genes.Expression profiling and enzyme activity of these families reveal their developmental and environmental influences.This study will pave the foundation to explore further the physiological roles of these members in MG detoxification as well as stress modulation of plants.

Identification and nomenclature of putative glyoxalase and D-lactate dehydrogenase members in tomato
To identify MG detoxifying members in tomato, previously characterized Arabidopsis glyoxalases (e.g.AT1G07645, AT1G06130, AT3G02720 etc.) and D-lactate dehydrogenase (e.g.AT4G36400, AT4G36400 etc.) proteins [9,17,26] were used as queries for Blastp search (e value: -1) against tomato genome 'Solanum lycopersicum ITAG3.2' in Phytozome v13 database [27].The resulting protein sequences were then analyzed using Pfam (e value: 1) [28] to confirm the presence of HMM profile PF00903, PF00753, PF01965.19 and PF01565 for GLYI, GLYII, GLYIII, and D-LDH, respectively.InterProScan [29,30] and SMART [31,32] were used to recheck the existence of the essential domains.The identified members were designated based on the sequence similarity with Arabidopsis orthologues by following the Arabic number system in alphabetic order and denoted by the prefix "Sl" for Solanum lycopersicum.Protein features of the identified members were predicted using the expassy protparam tools [33].The Plant-PLoc server [34] was utilized to determine the subcellular localization of the proteins.

Prediction of chromosomal localization and gene duplication events
Using information from the Phytozome Database, the locations of all newly identified glyoxalases and D-lactate dehydrogenases were mapped onto the twelve chromosomes of the tomato.Duplication occurrences were predicted using Blastp search (e-value 10 −10 ) in phytozome with 80% sequence identity.Tandemly duplicated (TD) genes were classified as those that included more than one homologous gene inside a 100 kilobase pair region on the same chromosome, whereas segmental duplicated (SD) genes were characterized as those that were not contained within the 100-kb region.The rate of synonymous substitution (dS) and non-synonymous substitution (dN), as well as evolutionary pressure (dN/dS) between the duplicated glyoxalase and D lactate dehydrogenase gene pairs, were determined by the PAL2NAL tool [35].Each duplicated gene pair's divergence period (T) was measured by T = dS/2λ, with a constant rate of 1.5108 substitutions per site per year for dicotyledonous plants.

Multiple sequence alignment and phylogenetic analysis of MG detoxifying members
Complete amino acid sequences of the putative glyoxalases and D-lactate dehydrogenases were aligned with their known orthologues from diverse plant species using ClustalW to identify the evolutionary relationship.Jalview 2 [36] was used to visualize the alignment.The evolutionary relationship and divergence of proteins were observed through the phylogenetic tree, built by MEGA-X [37], following the Neighbourhood-joining method with 1000 bootstraps.iTOL [38] was used for the modification and visualization of the tree.

Genomic structure analysis and protein architecture of the identified members
The exon-intron structures of glyoxalase and D-LDH genes were analyzed by comparing genomic and coding DNA sequences (CDS) through the Gene Structure Display Server (GSDS 2.0) [39].Gene features such as UTR, exon-intron, and intron phase locations were extracted from the GSDS server.Pfam and SMART were utilized to identify the domain positions and extract the amino acid sequences for the respective domains.Schematic representations of domain organization were manually created, and tables summarizing the presence of conserved residues within these domains were compiled.

Analysis of expression pattern of putative glyoxalases and D-LDH genes
Microarray data for tissue-specific expression of glyoxalases and D lactate dehydrogenase transcripts during various developmental stages of tomato were obtained through the Genevestigator Affymetrix tomato genome array database [42].Normalized Expression profile was obtained for seventeen different tissue types, including underground, aerial and reproductive tissue, as well as six development stages of vegetative growth and maturation-"shoot growth, inflorescence, flowering, fruit formation, ripening and complete ripening".Expression patterns in response to different biotic (pathogenic infection) and abiotic stressors (drought, salinity, heat) were retrieved from the same source.Genevestigator was used to calculate the log2 fold change ratio of each transcript under stress conditions as compared to the respective control.Heatmap was generated from the normalized expression data using the MeV software package [43] with the Manhattan distance correlation metrics [44].

Analysis of cis-acting regulatory element in the promoter regions
The putative promoter sequence (1 kb upstream of the start codons) of tomato glyoxalase and D-lactate dehydrogenase genes was obtained from the Phytozome database and analyzed through PlantCARE web tools [45] to identify the presence and type of cis-acting regulatory elements.The results of the analysis were represented in a bar diagram.The function of cis-acting elements was determined through a literature review.

Plant material and stress treatment
One of the Bangladeshi tomato varieties (BARI-4), seeds was collected from the Bangladesh Agricultural Research Institute (BARI), Bangladesh.Tomato seedlings were cultivated in a control condition with 14 hours of light and 8 hours of darkness at a temperature of 26±2˚C.The 8-day-old seedlings were given a variety of treatments.150 mM NaCl solution, 100 mM mannitol solution, and 5 mM H 2 O 2 solution were used to mimic salinity, drought, and oxidative stress, respectively, while plants were kept in 4˚C water or in an incubator at 40˚C to impose the cold or heat stress, respectively.Plants were kept under normal conditions for control.Shoot tissues were harvested at 0 h, 6 h, 12 h, and 24 h after the stress was applied.

Extraction of RNA extraction, and quantitative real-time PCR analysis
Total RNA was extracted using the 6 h post-treatment harvested frozen shoot samples using the manufacturer's instruction (SV Total RNA Isolation System, Promega Corporation, USA).The isolated RNA was quantified using a Thermo Scientific NanoDrop instrument and the first-strand cDNA was synthesized by using 10 μg of RNA according to the GoScript™ Reverse Transcriptase protocol (Promega Corporation, USA).Selected gene-specific primers were designed using the Primer-BLAST (https://ncbi.nlm.nih.gov/tools/primer-blast) and synthesized from Macrogen (https://dna.macrogen.com)along with SlEF1α as a reference gene [46] (S1 Table ).Real-time PCR was performed in the QUANTSTUDIO 1 3 REAL-TIME PCR SYSTEM with SYBR Green qPCR kits, and the fold change in expression for each gene at each condition was calculated using the 2 -ΔΔCt method [47].The experiments were carried out thrice for each condition and treatment.

Enzyme activity of GLYI and GLYII in response to different abiotic stresses
The total plant protein was extracted using the extraction buffer containing 0.1 M sodium phosphate buffer, pH 7.0, 50% glycerol, 16 mM MgSO 4 , 0.2 mM PMSF and 0.2% polyvinylpolypyrrolidone [19]; and quantified using Bradford method [48].The activity of GLYI was measured in 0.1 M sodium phosphate buffer (pH 7.5) with 3.5 mm MG and 1.7 mm reduced glutathione (GSH) as substrate by taking readings at 240 nm for the formation of S-Lactoylglutathione (SLG) over 5 min [19].Similarly, GLYII enzyme activity was measured by monitoring the reduction of SLG in 10 mm MOPS and 300 mM SLG [19].The activity of both enzymes was presented as μmol/min/mg of total protein.All the kinetic experiments were performed three times and the data were represented as the mean ± standard deviation (n = 3).

Ethics approval
All experimental research on plants, including the collection of plant materials has been used by following the relevant institutional, national, and international guidelines and legislation.

Identification of glyoxalases and D-lactate dehydrogenase members in tomato
Genome-wide analysis in tomato identified 13 putative GLYI, 4 GLYII, 3 unique GLYIII (DJ-1) and 4 D-LDH genes.The HMM Profile search using the Pfam database (E value 1.0) for the glyoxalase domain led to the pinpointing of 13 potential SlGLYI proteins with different transcript lengths (Table 1).There was no alternatively spliced transcript.They were designated as SlGLYI-1 to SlGLYI-11 based on their chromosomal position and similarities with Arabidopsis counterparts.Similarly, a total of 4 GLYII genes and 3 unique DJ-1 genes were identified by profile HMM search.GLYII genes were assigned as SlGLYII-1 to SlGLYII-3 based on similarity score and chromosomal position, while GLYIII (DJ-1) genes were labelled as SlDJ-1A to SlDJ-1D by following the previous report [9].
Initial screening for D-LDH found 24 genes having FAD_binding _4 domains.While these proteins had a common FAD-binding region, their catalytic activity might differ due to distinct secondary domains.Therefore, the domain organization of all identified FAD-binding superfamily members was analyzed to identify D-LDH, which specifically has D-lactate catalytic activity.In a previous study on sorghum, proteins containing only FAD_binding _4 domain as well as an additional FAD-oxidase_C were recognized as D-lactate catalyzing D-LDH, based on multiple sequence alignment and phylogenetic analysis [49].Thus, four D-LDH genes were predicted in our investigation encoding functionally active D-LDH, with two of them containing additional FAD-oxidase_C domain.They were numbered as SlDLDH-1-SlDLDH-4 (Table 1) according to blast hits against Arabidopsis thaliana D-LDH genes.

A comprehensive analysis of the identified SlGLYI and SlGLYII, SlDJ-1 and SlDLDH family members
The identified glyoxalases and D-LDH members were examined further to ascertain their physiochemical characteristics.The CDS length of SlGLYI members ranged from 414 bp (SlGLYI-  1).The isoelectric points (pI) of the SlGLYI family ranged from 4.92 (SlGLYI-4B) to 7.55 (SlGLYI-7 and SlGLYI-6A), indicating that cationic and anionic SlGLYI proteins could coexist simultaneously under certain physiological conditions.The bulk of SlGLYI proteins were localized in the chloroplast, with a few members detected in the plasma membrane, nucleus, mitochondria, and cytosol (Table 1).
Similarly, the coding DNA sequence of SlGLYII ranged from 723 bp (SlGLYII-1B) to 1461 bp (SlGLYII-3A), with an average length of 1040 bp.SlGLYII-1B was the largest member in terms of polypeptide length (486 aa) and molecular weight (54.5 kDa), while SlGLYII-3A was the smallest.Both the SIGLYI and SIGLYII protein families had positively and negatively charged proteins.Most SlGLYII proteins, like SlGLYI members, showed chloroplast localization except for SlGLYII-3B, which was found in the cytoplasm.
The physicochemical properties of all the identified GLYIII-like proteins (DJ-1) were also investigated as with traditional glyoxalase enzymes.The CDS length of SlDJ-1 ranged from 1167 (SlDJ-1D) to 1452 bp (SlDJ-1A).The smallest SlDJ-1 member, SlDJ-1D, had a length of 388 aa with a weight of 42.02 kDa, while the largest, SlDJ-1A, was comprised of 483 aa residues with a molecular weight of 51.74 kDa.The pI value ranges from 5.48 to 8.44.The Ploc server predicted that all the identified members were located in the chloroplast.
Similarly, D-LDH had CDS lengths ranging from 1701 bp (SlDLDH-4) to 1878 bp (SlDLDH-3).The identified genes have a mean molecular weight of 66 kDa.The pI of the majority of the identified members is acidic, except SlDLDH-3, with a basic pI.Likewise, most of the proteins were found to be localized in the cytoplasm, except SlDLDH-1 in the plasma membrane.

Chromosomal localization and gene duplication analysis
Glyoxalases and D-LDH genes were randomly dispersed throughout nine chromosomes of the tomato, whereas no genes were detected in the other three chromosomes (Fig 1B).On chromosome 2, there were a maximum of five methylglyoxal detoxifying members, followed by four each in chromosomes 1 and 12; three genes each in chromosomes 6 and 7; two genes in chromosome 11, and one gene each in chromosome 3, 8 and 9. Three SlDJ-1 genes were found to be distributed irregularly on chromosomes 1, 2, and 7.The D-lactate dehydrogenase genes were split across two chromosomes, with three genes on chromosome 2 and one on chromosome 11 (Fig 1B).
Three duplication events were observed among the identified genes.Between SlDLDH-3 and SlDLDH-4, there is a single duplication event on chromosome 2.A duplication between SlGLI-7A and SlGLYI-4A linked chromosomes 3 and 12, while another duplication between SlGLYII-3A and SlGLYII-3B linked chromosomes 12 and 7 (Fig 1B).To study the selection pressure on the duplicated SlGLYs and SlDLDH gene pairs, the value of non-synonymous (dN) and synonymous (dS) substitutions were estimated to discover whether the ratio is greater, less or equal to 1, indicating a positive, purifying and neutral selection respectively [50][51][52].All of the discovered duplicated gene pairs had a dN/dS ratio below 1, suggesting that purifying selection may play a role in their evolution (S2 Table ).Additionally, the divergence period for the duplicated gene pairs is estimated to be between 23.6 and 37.9 million years ago.

Gene and protein structure of conventional glyoxalases
The Gene Structure Display Server (GSDS 2.0) was used to depict the exon-intron structure of the genes to get an insight into the association between gene structure and activity.The exonintron patterns of the SlGLYI genes showed a high degree of variation.SlGLYI-3, SlGLYI-6B, and SlGLYI-6A all contained eight exons and distinct introns, but other members have different types of intron-exon organization (S1A Fig) .Genes encoding for GLYI enzymes fell into two groups based on sequence homology that exhibited distinct structural properties within each cluster (S1A Fig).By contrast, SlGLYII genes did not have such distinctive exon-intron combinations.SlGLYII-3 and SlGLYII-4 genes possessed nine exons and eight introns, whereas SlGLYII-3A contained nine exons and eight.SlGLYII-2 had 12 exons, the most of any SlGLYII gene, in both spliced versions (S1C Fig).
Domain architecture of putative SlGLYI elucidates that only the glyoxalase domain was present in all identified SlGLYI proteins (S1B Fig) .Three GLYI proteins, notably SlGLYI-3, SlGLYI-6A, and SlGLYI-6B, included two PF00903 domains, whereas the remaining proteins contained only one.The variation in the length of the glyoxalase domain has been linked with metal dependency [53,54].Thus, domain organization patterns may be used to predict the metal ion dependencies of SlGLYI proteins (S1B Fig) .Moreover, while all SlGLYII proteins had metallo-lactamase domains (S1D Fig) , only SlGLYII-1A possessed an additional HAGH_C (PF01623) domain, which was found responsible for the hydrolysis D-lactoyl-glutathione to form glutathione and D-lactate in human [55].Also, crystal structure analysis revealed that substrate binding might occur at the interface between HAGH_C and the catalytic β-lactamase region [56].Further, all identified SlGLYII contained the active site motif (C/ GHT) in addition to metallo-lactamase and thus was functionally active (S2B Fig).

Phylogenetic relationship, multiple sequence alignment and functional prediction of conventional glyoxalases
Phylogenetic analysis was used to examine the sequence similarities and evolutionary divergence of GLYI and GLYII enzymes of tomato along with various other reported species, including Arabidopsis, rice, G. max, Medicago, Sorghum, and Brassica.The GLYI proteins were found to be separated and clustered into three major clades (Fig 2A).Clades I and II represented GLYI enzymes with Ni 2+ and Zn 2+ dependency, respectively.In contrast, clade III is formed by inactive/functional diverge GLYIs.Clade-I comprised of Ni 2+ -dependent rice GLYI proteins, including OsGLYI-2, OsGLYI-7, and OsGLYI-11 and Arabidopsis AtGLYI-3 and AtGLYI-6 [53,57], suggesting Ni 2+ -dependent GLYI activity for members of SlGLYI in clade-I.Similarly, OsGLYI-8 and AtGLYI-2 in clade-II had Zn 2+ -dependent GLYI activity, indicating that other members of clade-II could be Zn 2+ -dependent enzymes [54].Therefore, based on the members in these clades, the Ni 2+ -dependent GLYI isoform may be more prevalent in plants than the Zn 2+ -dependent isoform (Fig 2A ).
Catalytic activity and metal dependency could be inferred for SlGLYI proteins by aligning their glyoxalase domain (PF00903) with Ni 2+ -dependent OsGLYI-11.2and Zn 2+ -dependent OsGLYI-8 orthologues [53,57], identifying the conserved actives site residues (S2A Fig) .Conserved residues H/QEH/QE, present in both OsGLYI-8 and OsGLYI-11.2,defined the metalbinding site of GLYI enzymes, where glutamate residues receive proton from the substrate [53,58].Only four members, namely SlGLYI-2, SlGLYI-3, SlGLYI-6A, and SlGLYI-6B possessed all the conserved residues in their N terminal glyoxalase domain and could predicted as active GLYI members (Table 2).According to previous studies, proteins with two GLYI domains of approximately 120 amino acids each were the putative Ni 2+ -dependent forms, whereas those with a single GLYI domain of approximately 142 amino acids and two additional sequence stretches were considered the putative Zn 2+ -dependent [53,57].Based on this, SlGLYI-2 could be Zn 2+ dependent, while the other three could be regarded as Ni 2+ dependent.

Phylogenetic relationship, multiple sequence alignment and structural features of tomato D-LDH
The exon-intron patterns in identified SlDLDH genes were entirely distinct from one another (Fig 4A).SlDLDH-1 possessed the greatest number of introns and exons, whereas SlDLDH-2   possessed the fewest.Both of these proteins had a FAD-oxidase _C domain and a FAD_bind-ing_4 domain.On the other hand, the remaining SlDLDH-3 and SlDLDH-4 proteins comprised the FAD_binding_4 domain solely (Fig 4B).Furthermore, phylogenetic analysis revealed that SlDLDH-1 and SlDLDH-2 cluster with AtDLDH, indicating that they are functionally comparable and likely to be mitochondrial proteins (Fig 4C).SbDLDH-3 and SlDLDH-4 proteins were more comparable to rice OsDLDH in terms of sequence similarities.All probable D-LDH proteins from tomato were analyzed with that of rice and Arabidopsis by sequence alignment (Fig 4D).The existence of highly conserved catalytic domains, denoted by black stars, suggested that the putative proteins are active (S4 Table ).

Expression pattern analysis of glyoxalases and D-LDH in different developmental and anatomical stages
Putative gene expression patterns of SlGLYI, SlGLYII, SlDJ-1, and D-LDH genes were obtained from a publically accessible microarray database using Genevestigator software to understand better the role of glyoxalases and D-LDH at distinct developmental stages and anatomical tissues (Fig 6).Expression analysis indicated that different tissues continuously expressed SlGLYI genes at varying levels.However, Ni 2+ -dependent SlGLYI-6A's expression level was greater from seedling to flowering and shoot development stages (Fig 6A ), opposite of inflorescence stages.By contrast, the expression of putative Zn 2+ -dependent SlGLYI-2 was significantly increased in some developmental phases, except ripening and finished ripening.While all SlDJ-1 expressed abundantly at the ripening stage, all SlD-LDH expressed oppositely.The expression levels of SlGLYI-7A, SlGLYI-7C and Zn 2+ -dependent SlGLYI-2 were higher in underground tissue such as root and root tips.In contrast, Ni 2+ -dependent SlGLYI-6A and SlGLYI-6B showed greater expression in aerial tissue such as seedlings, shoots, leaves and lamina and cotyledon of the seed, while SlGLYI-2 showed a similar expression pattern in these tissues (Fig 6B).Among SlGLYII genes, only SlGLYII-3A was expressed greatly in root and root tips, while SlGLYII-1B and SlGLYII-1A showed elevated expression in pericarp and exocarp tissue of fruits, which was similar to the case of SlDJ-1D.Both SlDLDH-3 and SlDLDH-4 were highly expressed in underground tissue like root and root tips, aerial tissue like shoot, steams and leaves as well as pistil and carpel tissue of flower.

Expression profiling under different abiotic stress and pathogenic infections
The expression of all the identified SlGLYI, SlGLYII, SlDJ-1 and SlDLDH transcripts were also analyzed under various pathogenic infections and abiotic stresses to understand their function in stress modulation of tomato through MG detoxification.The transcripts were expressed differentially in response to various abiotic stresses (Fig 7A).SlGLYI-9, SlGLYI-2, and SlGLYI-7A were upregulated in response to drought stress but were downregulated when salinity was applied.SlGLYI-9 expression was found to decrease in wounded fruits in either mature or ripening stages and when subjected to ammonium stress.Similarly, SlGLYII-3A and SlGLYII-B were upregulated during drought stress, and the expression somewhat decreased when exposed to salinity (Fig 7A).Among the analyzed SIGLYI, only SIGLYI-4B and SIGLYI-6A expressed at a reduced level in response to heat, and all SIGLYII except SIGLYII-3A expressed at a similar pattern.When exposed to heat, SIDJ-1D and SIDJ-1A were found to be upregulated.In response to local high and very high concentrations of ammonium, the transcript abundance of SlGLYII-6B, SlGLYII-3A, and SlDLDH-4 was highly expressed as compared to the other analyzed genes.Overall, most of the glyoxalase system enzymes were expressed at the highest level in response to drought stress, followed by heat and salinity.Most glyoxalase genes were upregulated in drought stress, while D-LDH were mostly upregulated in salinity and ammonium stress.
When tomato seedlings were subjected to a variety of biotic stresses, the SlGLYI gene expression pattern altered significantly (Fig 7B).SlGLYI-7A and SlGLYI-7C were highly upregulated in response to pathogen B. cinerea.However, their expression decreased when exposed to other pathogens.SlGLYII enzymes may play a larger part in the defence system since they are mostly upregulated in response to various pathogens, especially SlGLYII-1A and SlGLYII-1B.Among SlDJ-1 genes, only SIDJ-1D was elevated in tomato seedlings following infection with C. michiganensis, C. coccodes, G. intraradices, P. infestans, R. solanacearum, and T. urticae strains.Several SlGLYI and SlDLDH genes were more highly expressed, including SlGLYI-9, SlGLYI-2, SlGLYI-4B, and SlDLDH-4, whilst others exhibited little or no transcript alteration.

Validation of selected SlGLY, SlDJ-1, and SlDLDH genes under various abiotic stresses
The abiotic stress-responsiveness expression of the predicted enzymatically functional glyoxalase and D-DLDH genes were further studied in tomato using quantitative RT-PCR analysis under various stress conditions, including salinity, drought, oxidative, heat, and cold stresses (Fig 8 and S5 Table).All the predicted Ni 2+ -dependent SlGLYI genes namely, SlGLYI-3, SlGLYI-6A, and SlGLYI-6B were found to be more stress-inducible as compared with the only Zn 2+ -dependent SlGLYI-2 to all the stress treatments (Fig 8A).The putative Zn 2+ -dependent SlGLYI-2 showed marginal fluctuation in expression in all the stresses except for oxidative stress.Interestingly, three Ni 2+ -dependent members SlGLYI-3, SlGLYI-6A, and SlGLYI-6B behaved differently where SlGLYI-3 showed slight upregulation in expression in all the conditions, SlGLYI-6A showed more than 4 times upregulation, and SlGLYI-6B exhibited more than 6 folds downregulation.
Similarly, four functionally active GLYII genes namely SlGLYII-1A, SlGLYII-1B, SlGLYII-3A, and SlGLYI-3B also showed stress-mediated perturbations differently.Among them, two members such as SlGLYII-1A, and SlGLYI-3B showed more than 2 folds in expression in all five stress conditions, while the other two members SlGLYII-1B and SlGLYII-3A showed more than 6 folds of downregulation in expression (Fig 8B).Similarly, the expression of three noble GLYIII-like genes SlDJ-1A, SlDJ-1B, and SlDJ-1C was analyzed in tomato.SlDJ-1A showed significant upregulation under all these stress conditions, while SlDJ-1C showed significant downregulation (Fig 8C).The other member, SlDJ-1B showed a mixed pattern of expression where salinity, oxidative and cold stresses brought slight upregulation and drought and heat stresses induced downregulation of the gene (Fig 8C ).
Finally, expression profiling of all four identified SlDLDH isoforms was performed under five abiotic stress conditions (Fig 8D).The alteration in the expression of SlDLDH-1 and SlDLDH-3 was marginal as compared to the other two genes.SlDLDH-2 showed downregulation in the expression of more than 4 folds, while the expression of SlDLDH-4 was more than 4 folds upregulated under all stress conditions, potentially reflecting its role in metabolizing Dlactate during stress conditions.

Identification of putative cis-regulatory elements in the promoter region
Cis-regulatory elements were identified in the promoter region of the glyoxalases and D-LDH genes to understand their role in transcriptional regulation during developmental and environmental stimuli (Fig 9).The 1000 bp upstream region of SlGLY and SlDLDH transcripts were retrieved from the phytozome database and scanned to identify cis-elements in the promoter.Various hormone-responsive elements in the identified genes included abscisic acidresponsive (ABRE), ethylene-responsive (ERE), gibberellin-responsive (GARE), salicylic acidresponsive elements (TCA-element), and MeJA-responsive element (TGACG-motif) were  ).The positive alteration of the SlGLY and SlDLDH expression and their enzymatic activity under various abiotic stress conditions could be directly correlated with the presence of such a broad array of hormones and stress-inducible cis-elements.

Effect of various abiotic stress on the total activity of GLYI and GLYII enzymes
To gain deep insights into the function of glyoxalase genes in the abiotic stress adaptation in tomato, the total GLYI and GLYII activity was measured in response to various abiotic stress conditions including salinity, osmotic stress, heat, cold, heat, and drought, and compared with the respective untreated control conditions (Fig 10 and S7 Table ).The plants were grown normally and exposed to different stresses in the laboratory.The activity was measured in four different time points such as 0h, 6h, 12h and 24h and expressed as μmol/min/mg of total protein.An enhancement of total GLYI activity was observed for salinity and osmotic stress.However, the activity decreased over time (Fig 10A).On the other hand, heat, cold and drought negatively affected the total GLYI activity, which gradually decreased over time.A similar pattern was observed for the total GLYII activity where salinity, osmotic and heat (to some extent) enhanced the activity (Fig 10B).Both the activity of GLYI and GLYII enzymes were found to

Discussion
Climate change, which results in various forms of abiotic stresses, is the main contributor to yield losses in various economically significant crops worldwide.Crops prone to dryness, high temperatures, extreme cold, and saline soil as a result of global warming suffer considerable yield losses.The production of tomato has continuously declined in the last century due to adverse climatic conditions.The metabolic byproduct MG is produced spontaneously in all living organisms [6], while the level goes up in response to diverse abiotic stressors [20].The importance of the glyoxalase system in eliminating MG and providing tolerance from a variety of abiotic stressors has already been established through previous research [5,6].Previous studies had also identified glyoxalase genes in different plant species, including Arabidopsis (10 GLYI, 5 GLYII, and 6 DJ-1) (8,26), Soybean (24 GLYI, 12 GLYII, and 7 DJ-1) [14], M. truncatula (29 GLYI, 14 GLYII, and 5 DJ-1) [60] and rice (11 GLYI, 3 GLYII, and 6 DJ-1) [26].Four D-LDH genes were reported in Sorghum bicolor, with 15 GLYI and 6 GLYII genes [49].This study identified 13 putative GLYI, 4 GLYII, 3 DJ-1, and 4 DLDH genes in the tomato genome Solanum lycopersicum ITAG3.2 (Table 1).
Among the identified members in tomato, three duplication events were observed in this study (Fig 1B).Two duplication events were discovered among the conventional glyoxalase genes (GLYI and GLYII), whereas one duplication was found in the SlDLDH genes.All duplication events were determined to be segmental (S2 Table ).The SlGLYI-4A/7A pair was formed 37.9 Mya ago, whereas SlGLYII-3A/3B and SlDLDH-3/4 were formed 23.6 and 29.5 Mya ago, respectively.In contrast, soybean exhibits ten duplication patterns among glyoxalase genes (GLYI and GLYII) [14], whereas five GLYI and one GLYII gene have reported such occurrences in Medicago [60].These duplication events in plants have led to gene gain/loss and confer abiotic stress modulation [63].
The predicted SlGLYI proteins may be classified into two categories based on their activity, similar to their orthologues in other plant species.Active GLYI is either Zi 2+ -or Ni 2+ -dependent.The metal selectivity of SlGLYI proteins was estimated based on their domain sequencing and length, as indicated in the previous studies [57].Only four GLYI proteins, SlGLYI-2, SlGLYI-3, SlGLYI-6A and SlGLYI-6B, were considered active as they had conserved metalbinding sites (Table 2).SlGLYI-2 was discovered as a Zn 2+ -dependent glyoxalase protein and had a domain length of 143 amino acids (S1B Fig), whereas the remaining three proteins are likely to be Ni 2+ -dependent with domain lengths of approximately 120 amino acids, quite similar to its orthologues in rice [53,57].Other identified SlGLYI might be functionally inactive GLY-like proteins that lacked the conserved motif.
SlGLYII-3A and SlGLYII-3B were discovered as functionally active GLYII enzymes due to their conserved metal-binding site THHHXDH (Table 3) and high sequence similarity to rice's functionally active OsGLYI-2 and OsGLYI-3 proteins [20].SlGLYII-1A and SlGLYII-1B, on the other hand, were potential SDO enzymes based on sequence similarity to the OsGLYII-2 and AtGLYII-2 proteins [59] that exhibited SDO activity, indicating functional divergence among GLYII enzymes (S2B Fig) .In previous research, the Hydroxyacylglutathione hydrolase (HAGH_C) domain was also discovered in the C terminal of functionally active GLYII enzymes by sequencing and crystal structure analysis postulated that substrate binding might occur at the interface between HAGH_C and the catalytic β-lactamase region [56].As a result, the original GLYII enzyme may have an additional HAGH_C domain, as confirmed with SlGLYII-1A (S2B Fig) .In-silico investigation of SlDJ-1 proteins found three candidates in tomato with DJ-1 domain (Table 1), whereas five members were identified in Medicago [60], six in Arabidopsis and rice [8,26] and seven in soybean [14].Previous research on DJ-1 protein revealed plants had three major types of DJ-1 proteins, diverged from prokaryotes, protists and green algae [64].Most of the plant-derived DJ-1 proteins had double DJ-1/PfpI domain, except for date palm, which included up to four [65].In the glutathione-independent pathway, the DJ-1/PfpI containing glyoxalase III enzyme converts MG to harmless D-lactate in a single-step reaction [9].D-LDH enzymes belonging to the FAD_binding_4 superfamily catalyze D-lactate's transformation to pyruvate.In the HMM profile search, 24 such members were initially identified as belonging to the FAD-binding superfamily.In addition to FAD binding 4, these identified members featured several additional domains, demonstrating the functional variety of the group.To anticipate SlDLDH members that specifically catalyze D-Lactate transformation, SlDLDH containing a single FAD_binding_4 or an additional FAD-oxidase_C domain were selected according to a prior study [49].Therefore, four such enzymes with specific D-lactate catalytic activity were found in tomato, two of which had an additional FAD-oxidase C domain.Homology modelling of SlDLDH-1 revealed that it contained all the active site residues at R428, I446, H479, H486, E523, and H524 to catalyze D-lactate (Fig 5E ), indicative of functional activity.
Previous research on rice and Arabidopsis identified glyoxalase genes to be highly selective to certain tissues or developmental stages [26].In a similar study on M. truncatula, MtGLYI-4 showed the greatest expression level across all 17 analyzed tissues, while MtGLYI-3, MtGLYI-18, and MtGLYI-20 had the lowest expression levels [60].In the current study, we used publicly available microarray data and observed a cluster of GLYI genes, including SlGLYI-11, SlGLYI-6B, and SlGLYI-6A, expressed in great abundance in aerial tissue like shoot, leaf and lamina while showing minimal expression in underground tissues such as root, root tips (Fig 6B).Recent research has shown that one of the rice glyoxalases, OsGLYI-3, promotes seed lifespan and salt tolerance in rice plants [66].SlGLYI-11 may serve a similar function as it showed significant expression in seed tissue and was highly elevated under saline conditions.However, further experimental evidence for functional characterization.
Tomato GLYII member SlGLYII-3A was significantly expressed in root and root tips while SlGLYII-1B and SlGLYII-1Ashowed the highest expression level in fruits' Exocarp and Pericarp walls (Fig 6B).This demonstrates the presence of functional diversity among the glyoxalase II protein in various developmental cues; a similar phenomenon was observed in the case of soybean, where GmGLYII-12, GmGLYII-5 and GmGLYII-4 were found to be only expressed in seed tissue, while GmGLYII-9 was expressed in seeds, roots and flower tissue [14].SlDJ-1D showed the highest expression in fruit tissue while the expression of SlDJ-1C and SlDJ-1A was very minimal (Fig 6B).Similar phenomenon was observed in M. truncatula where MtDJ-1D and MtDJ-1A showed expression equally in underground, aerial and seed tissue, while MtDJ-1C and MtDJ-1B showed no expression [60].D-lactate dehydrogenase in Sorghum bicolor, SbDLDH-3 and SbDLDH-4 showed a higher level of expression in root and shoot tissue [49], similar to SlDLDH-3 in tomato, which is also highly expressed in the root (Fig 6B ), indicating its common role of root development in different plant species.The increased expression of SlGLYII-1B, SlGLYII-1A, and SlDJ-1D genes during the fruit ripening stage revealed that GLYII and DJ-1 m play may play some roles in fruit maturation of tomato (Fig 6A).However, the majority of GLYI and D-LDH may play a part in shoot formation, Inflorescence and flower development, as can be observed from the heatmaps (Fig 6A).
In previous studies on rice [26], G. max [14] and Medicago [60], different members of the glyoxalase family had been identified as potential modulators of abiotic stresses.In this study, the transcriptomic data from Genevestigator revealed that most of the glyoxalase genes were upregulated in response to drought, followed by heat and salinity.In most scenarios, upregulated genes are clustered together, a pattern that is also common for downregulated genes (Fig 7A).Members of the SlGLYI gene family-SlGLYI-9, SlGLYI-7A, along with SlGLYII-3A and SlDJ-1D were elevated in response to dehydration and heat stress but decreased in response to salinity exposure.The Ni 2+ -dependent SlGLYI-6A exhibited increased expression under salinity and drought conditions, whereas SlGLYI-6B was downregulated.Conversely, the Zn 2 + -dependent SlGLY-2 was primarily expressed in response to drought, heat, and physical stress (Fig 7A).
In previous research, the DJ-1 members of prokaryotes Escherichia coli and eukaryotes S. cerevisiae displayed salt tolerance [67], which may be confirmed by the upregulation of tomato SlDJ-1A and SlDJ-1C in response to salinity (Fig 8C ).A similar phenomenon was observed in Medicago, in which all of the MtDJ-1 members were expressed under salinity stress along with other stress conditions [60].A recent investigation into tomato glyoxalase III (GLYIII) genes revealed that their overexpression enhances salt and osmotic stress tolerance, contributing positively to plant growth and yield [68]. Further, transgenic tobacco expressing E. coli DJ-1 demonstrated enhanced salt tolerance during vegetative and reproductive stress.Transgenic sugarcane overexpressing GLYIII (EaGlyIII) also improves salt and drought tolerance, according to a recent study [69].Among SlDLDH members, SlDLDH-4 was mostly upregulated under abiotic stress conditions while SlDLDH-2 was highly downregulated.In contrast SlDLDH-1 and SlDLDH-3 showed heterogenous expression in response to abiotic stress (Fig 8D ).
Moreover, microarray data revealed transcript abundance of SlGLYI-7A and SlGLYI-7C of the GLYI family, SlGLYII-1B and SlGLYII-1A of the GLYII family, along with SlDJ-1C and SlDLDH-4 against the pathogenic fungal infection of Botrytis cinerea.Furthermore, most of the SlGLYIIs with SlGLYI-7A, SlGLYI-2, SlDJ-1D, and SlDLDH-4 showed greater responses against various pathogenic infections (Fig 7B).Alongside with qPCR data, the activity of total GLYI and GLYII enzymes was found to gradually decrease with the increment of stress duration, indicating the physiological response of tomato glyoxalase in stress modulation ( Fig 10).Our study on the functional characterization of MG detoxifying genes could lead to the development of a stress-tolerant tomato variety with improved pathogen resistance.

Conclusion
To summarise, genome-wide investigation of methylglyoxal detoxifying enzymes in tomato led us to identify 13 GLYI, 4 GLYII, 3 DJ-1, and 4 D-LDH genes distributed in 9 different chromosomes.Three segmental duplication events were observed.Based on the previous literature study, sequence analysis, the presence of conserved domain architecture, the metal ion dependency and putative enzymatic activity of the identified proteins have been predicted.The expression data confirmed the significant role of these genes in different developmental stages, anatomical tissues, and stress modulation pathways by mitigating MG accumulation during different biotic and abiotic stresses.SlGLYI-6A, SlGLYII-3B, SlDJ-1A and SlDLDH-4 of the MG detoxifying gene family showed significant upregulation in response to environmental stressors, while SlGLYII-1B, SlGLYII-1A, SlGLYI-2, SlDJ-1D, SlDLDH-4 were highly expressed in response pathogens, indicating their role in plant defense systems and stress modulation.Further functional analysis of these genes may help to develop an engineered stress-tolerant plant.

Fig 1 .
Fig 1. Role of methylglyoxal detoxifying members in plants and their chromosomal location in tomato.(A) Schematic diagram of the glyoxalase route for methylglyoxal (MG) detoxification in plants.In the glutathione-dependent pathway, MG is converted to S-D-lactotyl glutathione and D-lactate sequentially by GLYI and GLYII.GLYIII catalyzes the reaction in a single step without the aid of reduced glutathione.D-lactate is then converted to pyruvate by CYTc-dependent D lactate dehydrogenase (DLDH).(B) Chromosomal distribution of glyoxalases and DLDH genes in Solanum lycopersicum.The presence and position of all identified genes were shown in nine different chromosomes with three segmental duplication events added by a straight line.https://doi.org/10.1371/journal.pone.0304039.g001

Fig 2 .
Fig 2. Evolutionary relationship of conventional glyoxalase proteins.The phylogenetic tree was constructed for the (A) GLYI and (B) GLYII proteins from Sorghum, rice, Arabidopsis, Medicago, tomato, Brassica, soybean, and other plants using the Neighbour-Joining method in MEGA X with 1000 bootstrap replicates.Different clades indicated the different subgroups based on co-factor dependency and functional similarity.https://doi.org/10.1371/journal.pone.0304039.g002 Phylogenetic relationship, multiple sequence alignment and structural features of tomato DJ-1 membersGSDS 2.0 was used to analyze the gene features of all identified DJ-1.SlDJ-1A possessed the most CDS regions, roughly ten, with a single upstream/downstream sequence.About five CDS regions were identified in the SlDJ-1D gene, compared to SlDJ-1C, which had eight.SlDJ-1A was the longest gene, whereas SlDJ-1C had the smallest (Fig3A).In previous studies, DJ-1 proteins of plants such as Arabidopsis and rice were reported to have two DJ-1/PfpI domains, whereas humans, Drosophila, and E. coli only had a single domain[9].The DJ-1/PfpI domain contained roughly 140 to 150 amino acids in most organisms, with the only exception being E. coli.According to multiple sequence alignment, most SlDJ-1 proteins had two consecutive DJ-1domains (Fig3B) except for SlDJ-1A, which contained three DJ-1/PfpI domains of 70, 105, or 165 amino acids, respectively.

Fig 3 .
Fig 3. Phylogenetic relationship, sequence alignment, gene structure and domain architecture determination of SlDJ-1 members.(A) Gene structure was drawn by GSDS 2.0.(B) Schematic representation of domain organization of SlDJ-1 proteins.(C) A phylogenetic tree was built for the evolutionary conservancy among the DJ-1 members from various plant species by MEGA-X using the Maximum likelihood method with 1000 bootstrap.(D) Multiple sequence alignment was marked with black stars to indicate the conserved active site residues.https://doi.org/10.1371/journal.pone.0304039.g003

Fig 4 .
Fig 4. Elucidation of the structure and phylogenetic relationship DLDH members.(A) The exon-intron structure of SlDLDH genes.(B) Schematic representation of the domain architecture of SlDLDH proteins.(C) The evolutionary relationship of SlDLDH proteins with their orthologue from Sorghum, rice and Arabidopsis was observed by constructing a phylogenetic tree using MEGA-X with the Maximum Likelihood method and 1000 bootstrap.(D) Sequence conservance of SlDLDH proteins to identify the evolutionary conserved active residues that were marked with downward arrows.https://doi.org/10.1371/journal.pone.0304039.g004

Fig 6 .
Fig 6.Expression profiling of methylglyoxal detoxifying genes from tomato.Expression profiling of glyoxalases and DLDH genes of tomato was performed with hierarchical clustering at (A) different developmental stages and (B) different anatomical tissues.Expression data were obtained for 17 different tissues and 6 distinct developmental stages from the Genevestigator database.The MeV software package (http://mev.tm4.org/) was used to build a heatmap with hierarchical clustering based on the Manhattan correlation.The level of expression can be interpreted by the colour scale where the intensity of colour is directly proportional to the level of expression.https://doi.org/10.1371/journal.pone.0304039.g006

Fig 7 .
Fig 7. Expression analyses of tomato glyoxalases and DLDH genes in response to various abiotic and biotic stresses.Relative expression data of all the available SlGLYIs, SlGLYIIs, SlDJ-1s, and SlDLDHs were retrieved from the genevestigator and combined and analyzed in response to different types and duration of (A) abiotic and (B) biotic stresses.Multiple Experiment viewer software packages (http://mev.tm4.org/) were used to construct the heatmaps with hierarchical clustering by using the Manhattan distance metric.The yellow squares in the heatmap indicate up-regulation and green indicates down-regulation of the corresponding transcripts.https://doi.org/10.1371/journal.pone.0304039.g007

Fig 8 .
Fig 8. Abiotic stress-responsiveness of putative functional MG detoxifying members in tomato by qRT-PCR.The expression pattern of in silico analysed putative functional four SlGLYI (A), four SlGLYII (B), three SlDJ-1 (C), and four SlDLDH (D) genes was investigated in 10d old tomato seedlings by exposing different abiotic stress including salinity, drought, oxidative, heat, and cold stress for 6 h.Bar diagrams illustrate the average fold changes in gene expression relative to controls using SlEF1α as a reference gene.Each bar is distinctly coloured based on the type of stress conditions.https://doi.org/10.1371/journal.pone.0304039.g008

Fig 9 .
Fig 9.In silico promoter analysis of SlGLYI, SlGLYII, SlDJ-1 and SlDLDH members.The 1 kb upstream sequence of all the candidate genes was retrieved from the phytozome database and scanned through PlantCARE.The core promoter elements such as TATA-box and CAAT-box were not shown in the figure for better visibility.The number of stress-related and hormonal-related -cis-elements was determined for the analyzed promoter sequences.https://doi.org/10.1371/journal.pone.0304039.g009

Fig 10 .
Fig 10.Measurement of total GLYI and GLYII activity in response to various abiotic stresses.Total GLYI (A) and GLYII (B) enzyme activity were measured in response to various abiotic stresses such as salinity, osmotic, heat, cold, heat, and drought stresses at four different time points of stress exposure (0h to 24h).The activity was represented as μmoles/min/mg of total protein.All the experiments were repeated thrice and represented as the average ± standard deviation (n = 3).https://doi.org/10.1371/journal.pone.0304039.g010

Table 1 . List of identified glyoxalases and D-lactate dehydrogenase members in tomato along with their detailed information and localization. Sl no Gene Name
Subsequently, SlGLYI-6A was the largest protein of the SlGLYI family with a length of 336 aa and size of 44.43 kDa, while SlGLYI -1 was the smallest one with 173 aa in length and 44.43 KDa in size (Table Abbreviations: CDS = Coding DNA sequence, MW = Molecular weight, bp = base pair, aa = Amino Acid, KDa = Kilo Dalton, pI = isoelectric point https://doi.org/10.1371/journal.pone.0304039.t0011) to 1191 bp (SlGLYI-6A).

Table 2 . Information on domain organization of SlGLYI proteins for the prediction of enzymatic activity and metal ion dependency. GLYI Proteins GLYI domain Presence of conserved metal-binding site Expected enzyme activity Metal ion dependency
[14,59,60] to have sulfur dioxygenase (SDO) activity[14,59,60], indicating that members of the GLYII family functionally diverged.The remaining clade consisted of the Glyoxalase II-like proteins (Fig2B).All SlGLYII were regarded as functionally active as they had both active site (C/GHT) and metal ion binding site THHHXDH (Table3), similar to that of AtGLYII-2 and OsGLYII-2, denoted with black lines in the multiple sequence alignment (S2B Fig). https://doi.org/10.1371/journal.pone.0304039.t002werepreviously