Bioremediation of Chromium-Contaminated Groundwater Using Chromate Reductase from Pseudomonas putida: An In Silico Approach

Abstract

Chromium is a toxic heavy metal abundantly present in the environment, specifically in groundwater. The groundwater in Saudi Arabia was assessed for heavy metal presence; chromium was detected at a high concentration in Madinah. Many researchers have used various bioprocesses over the last few decades to mitigate Cr(VI) toxicity. The genus Pseudomonas member Pseudomonas putida is widely dispersed in the natural environment. P. putida is chromate-resistant and has a high chromate reduction rate. Bioremediation procedures can eradicate the most potently toxic metal, Cr(VI), in water, air, and soil. Chromate reductase (ChrR) is a bacterial enzyme from P. putida that can be utilized in bioremediation to remove chromate from the environment in a cost-effective and environmentally safe approach. To comprehend the role of ChrR in reducing Cr(VI) to Cr(III), a thorough sequence analysis was followed by constructing models for wild-type and mutants by applying several homology modeling techniques. The protein structure quality of the generated models was evaluated, and the best model was adopted for further optimization by employing an energy minimization technique. Molecular docking studies investigated the intra-molecular interaction between wild-type and mutant ChrR and Cr(VI). Our study is a novel method for determining the 3D structure and interaction of ChrR with Cr(VI) to convert it to a less hazardous form (III). Additionally, it provides stable mutants: Arg83Trp, Gly124Ile, and His127Trp, with a high binding affinity for Cr(VI), which can be considered for protein engineering to produce stable and efficacious enzymes to reduce Cr(VI) to a less toxic form.

1. Introduction

Saudi Arabia’s economy is one of the fastest-growing in the Gulf, with thriving oil and gas, agricultural, and pharmaceutical industries. This increased urbanization, population, and demand for goods and services. All of these affect the country’s natural resources, especially water [1]. In KSA, seawater desalinization meets the needs of people in or near coastal areas, so most of the population relies on non-renewable resources [2]. Metal contamination of the environment poses a severe risk to human health on a global scale. The quality of groundwater is influenced by several factors, including how water interacts with sediments and soils, the flow pattern, the types of rocks, and the geochemical conditions that are common, such as dissolution, redox state, leaching, precipitation, and ion exchange [3]. It is known that the toxicity of metals is induced by the possible creation of highly reactive species, such as reactive oxygen species [4] and certain other free radicals, which can damage DNA, deplete proteins, and oxidize lipids, among other impacts on cell structure [5]. Researchers found that hazardous metal ions contaminate KSA’s groundwater in several areas. Saudi Arabia’s groundwater contains Al, As, Ba, Cd, Co, Cr, Cu, Fe, Hg, Li, Mn, Ni, Pb, Se, V, and Zn. Madinah, in Saudi Arabia, was the focus of the investigation for its historical importance. A quality assessment of the groundwater was conducted by Maghraby et al. in the southern Madinah region of Saudi Arabia. The authors selected the study region as being adjacent to Wadi al Aqiq, located about 20.0 km south of Madinah, and 29 water samples were collected. As: 1–2, Cd: 1–4, Cu: 2–8, Fe: 1–37, Cr: 1–146, Ni: 1–18, Mn: 1–39, and Zn: 22–475 µg/L were the metal ions found in the water samples of Madinah. According to the authors, Pb, Mn, Ni, Cu, As, and Zn concentrations were within acceptable levels, whereas Cd and Cr exhibited higher quantities [6]. Bamousa and El-Maghraby later performed a quality study and identified the sources of contamination in Madinah, Saudi Arabia. In the regions of Uhud, Quba, Al Aqool, and the vicinity of the Prophet’s Holy Mosque, the researchers analyzed 32 groundwater samples [1]. The reported metal ion concentrations were as follows: Fe: 0.046–0.67, Pb: 0.0015–0.027, Mn: 0.011–0.48, Zn: 0.011–0.29, Ni: 0.001–0.18, Cu: 0.0023–0.0087, Cr: 0.011–0.11, Cd: 0.0014–0.083, Se: 0.0001–0.017, As: 0.0001–0.045, V: 0.002–0.044, Hg: 0.0001–0.0007, Al: 0.0001–0.11, and Co: 0.0001–0.0014 mg/L. According to the WHO, several samples from the southwest were found to contain high amounts of Al, As, Cd, Cr, Ni, and Pb, although the authors reported normal proportions of Cu, Co, Fe, Hg, Mn, Se, V, and Zn metal ions in the water. Due to its potential for mutagenesis effects, Cr is classified by the International Agency for Research on Cancer (IARC) as a human carcinogen (Group 1 category) [7]. Drinking water is the main source of chromate exposure [8]. In addition, the National Toxicology Program conducted a cancer bioassay using chromate in mice and rats using drinking water in 2007 [9].

Microorganisms are primitive and adaptable. They adapt by rearranging DNA and swapping genetic material needed for nutrient cycling, primary production, and pollutant catabolism [10]. Arsenic, chromium, copper, cadmium, lead, nickel, mercury, molybdenum, zinc, and vanadium are frequent metal contaminants. Metal pollution causes chronic and degenerative diseases [11]. There are seven different oxidation states of chromium (0–VI), with the metallic Cr(0), trivalent [12], and hexavalent Cr(VI) states being the most abundant in the environment and industries [12]. The possibility that bacteria use chromium (VI) as a terminal electron acceptor when oxidizing organic compounds has recently attracted attention. According to published research, different taxonomic groups of bacteria cultures can reduce chromium (VI) by releasing Cr (OH)3. Additionally, biotechnologies for galvanic wastewater treatment have been developed and implemented. It has been demonstrated that a biotechnology approach may reduce the chromium (VI) concentration in wastewater to the nominal limit of 0.3 mg/l. Compared to the conventional physico-chemical treatment, which is widely employed, this biological method is both less expensive and less harmful to the environment [13]. Microorganisms resistant to chromium include B. subtilis, B. cereus, P. putida, Ps. ambigua, Ps. Aeruginosa, Ps. fluorescens, A. eurydice, E. cloacae, E. coli, M. roseus, D. desulfuricans, and D. vulgaris. Pseudomonas is a bacterial genus with a greater biodegradation capacity [14]. Pseudomonas sp. is among the most common Gram-negative microorganisms. Pseudomonas sp. is present in the soil, groundwater, vegetation, and domestic environments such as hot tubs. Numerous chemical substances, such as fatty acids, pesticides, and aliphatic and aromatic hydrocarbons, can be broken down by some Pseudomonas species [15]. By combining genetic engineering and bioremediation, it is possible to manipulate the bacteria genome, improving their ability to detoxify toxic metals. Single-gene or operon-based genetic engineering, pathway creation, and sequence-based modifications to pre-existing genes are all examples of such methods [16]. They are, therefore, advantageous for application as bioremediation agents.

In silico studies can quickly and simply explore enzymes and their activity, as opposed to time-consuming, labor-intensive, and expensive experimental procedures [4]. Analyses were conducted on chromate reductase (ChrR) from P. putida to determine the presence of motifs and domain, its family, and evolutionary connections. In the absence of the P. putida ChrR structure, homology models made from closely related structural homologs would be the perfect way to start comprehending the binding interaction with Cr (VI). The mutant models were studied to identify the most potent residues that are effective for protein engineering and modifying it to a stable and efficient enzyme for bioremediation. The research sheds light on the ligand-binding mechanism and the structural features of the ligands that control the biological activity of this bacterium. These findings contribute to the selection of P. putida for bioremediation methods. Future research on structural and dynamic aspects will aid in developing stronger chromate reductase-related enzymes that might be used to reduce chromate pollution in the environment effectively.

2. Materials and Methods

The computational framework created for sequence annotation is shown by stages I, II, III, IV, and V. Stage I includes the following components: retrieval of ChrR sequence, its characterization, comparison, prediction of essential function, and virulence factor. In stage II, reference templates for model development are found, and the generated models are then evaluated using various bioinformatics techniques. In stage III, docking investigations were carried out to comprehend how chromate reductase (ChrR) and hexavalent chromium (Cr(VI)) interacted. Stage IV included the identification of mutant residues, development of mutant models, docking analyses of the mutant models, and Cr(VI). Stage V involves the detection of a closed link between the heavy metal Cr(VI) and the receptor.

2.1. Analysis of ChrR Sequence

The National Center for Biotechnology Information’s (https://www.ncbi.nlm.nih.gov/ (1 June 2022)) online data retrieval system was used to download ChrR protein sequences from P. putida in FASTA format [17]. The physiochemical characteristics of ChrR protein sequences, including their molecular weight, extinction coefficient (κ), instability index, isoelectric point (pI), aliphatic index, and grand average of hydropathicity, were measured in great detail using Expasy’s ProtParam (https://web.expasy.org/protparam/ (1 June 2022)) service (GRAVY). A key step in determining the virulent protein that could be the target of the development of antimicrobial drugs is the identification of the virulence factor in proteins. Here, it is necessary to have a non-virulent protein capable of oxidizing Cr(VI) to Cr(III). Using the powerful algorithm of Virulentpred (http://203.92.44.117/virulent/ (1 June 2022)), VICMPred (https://webs.iiitd.edu.in/raghava/vicmpred/ (1 June 2022)), and MP3 (http://metagenomics.iiserb.ac.in/mp3/index.php (1 June 2022)), the virulence of the protein sequences was predicted. The program Virulentpred accurately predicts bacterial virulence factors using support vector machine (SVM) techniques [18]. The web server VICMpred predicts the function of Gram-negative bacterial proteins [19]. The MP3 program predicts pathogenic proteins in metagenomic datasets [20].

The retrieved protein sequence was submitted to the identifying of ortholog sequences, family, and domains, and therefore their functions were evaluated. The tool eggNOG 5.0 (http://eggnog5.embl.de/ (2 June 2022)) was used to identify orthologs of ChrR in 4445 species of bacteria [21]. PSI-BLAST (https://blast.ncbi.nlm.nih.gov/ (2 June 2022)) [22] was applied to the Protein Data Bank (PDB) database (https://www.rcsb.org/ (2 June 2022)) [23] to identify structurally related proteins.

Various bioinformatics methods were utilized for the exact attribution of protein functions. The protein sequence families were found using Pfam (http://pfam-legacy.xfam.org/ (3 June 2022)) [24] and CATH (https://www.cathdb.info/ (3 June 2022)) [25]. CDD (https://www.ncbi.nlm.nih.gov/cdd/ (3 June 2022)) [26], CDART (https://www.ncbi.nlm.nih.gov/Structure/lexington/lexington.cgi (3 June 2022)) [27], and SMART (https://smart.embl.de/smart (3 June 2022)) [28] online tools were adopted for predicting domains. CDD supports multiple sequence alignments of domains and whole proteins to identify conserved domains. The Simple Modular Design Research Tool (SMART) investigates and annotates the architecture of protein domains. Proteins with comparable sequences are grouped together by CDART, which also assigns a score based on the architecture of the protein. Consensus was determined by analyzing the data obtained from the aforementioned tools, and the ProFunc (https://www.ebi.ac.uk/thornton-srv/databases/profunc/ (4 June 2022)) tool’s [29] outcomes were taken into account for functional assignment. Protein sequence motifs are thought to be the distinguishing characteristic of protein families, making it easier to predict how proteins may function. Since they are connected to catalytic processes, motifs play a significant role in enzyme function. InterProscan (http://www.ebi.ac.uk/interpro/ (4 June 2022)) [30] was used to identify motifs in the query sequence. InterProscan is comprised of several well-developed algorithms for protein signature recognition, such as PANTHER, ProDom, Pfam, PRINTS, Prosite, and SUPERFAMILY. MEME (https://meme-suite.org/meme/ (4 June 2022)), which employs an expectation-optimizing algorithm to fit a two-component finite hybrid algorithm to a set of sequences, and the Motif search tool (https://www.genome.jp/tools/motif/ (4 June 2022)) were also used for motif searching.

2.2. Multiple Sequence Alignment and Evolutionary Analysis

Multiple sequence alignment (MSA) serves a key role in clarifying sequence relationships and identifying functionally significant conserved areas [31]. MSA, ClustalW (https://www.genome.jp/tools-bin/clustalw (7 June 2022)) [32], and MEGA version 5.1 (https://www.megasoftware.net/ (7 June 2022)) were applied to identify conserved areas in the sequence [33]. The neighbor-joining algorithm (N-J algorithm) and Poisson substitution approach were employed in phylogenetic analysis with MEGA version 5.1 to deduce sequence homology and establish evolutionary relationships.

2.3. Secondary Structure Assessment

Secondary structure prediction is a crucial stage in the identification of protein characteristics and fold recognition, as well as an intermediate step in the prediction of three-dimensional structures. Since PSIPRED (http://bioinf.cs.ucl.ac.uk/psipred/ (8 June 2022)) employs neural network methods based on position-specific score matrices produced by PSI-BLAST, it was chosen as the preferred method for predicting secondary structures of the query sequence [34]. This approach of prediction is superior to all others, as it has been demonstrated to yield the highest Q3 score.

2.4. Three-Dimensional Structure Generation

A tertiary structure for the query sequence was built using an in silico approach with precision comparable to the experimental data in order to demonstrate the structure–function link. The homology modeling method, which constructs the structure based on a reference template with a known experimentally established tertiary structure, was used to develop the structures. Structures are more evolutionary conserved than sequences, and similar sequences adopt similar structures. Thus, the produced 3D structures will establish the groundwork for comprehending the novel functionalities of the protein sequences that have not yet been annotated and offer substantial insight into interaction investigations.

The PSI-BLAST program was used to carry out a sequence similarity search against the PDB database to provide a list of results that were most similar to the query sequence in order to identify the template. The sequence with the highest sequence identity, the largest query coverage, and the score with the lowest e-value was chosen as the reference template. Based on the outcomes returned by BLAST, the coordinates of the template structure were obtained from the Protein Data Bank [23].

2.5. Three-Dimensional Model Generation and Optimization

The five homology modeling tools listed below were selected to construct a three-dimensional structure: (i) Build Homology Model protocol from Biovia Discovery Studio v21.1 (Dassault Systèmes BIOVIA, Discovery Studio Modeling Environment, Release 21.1, San Diego: Dassault Systèmes, 2021), (ii) SwissModel (https://swissmodel.expasy.org/ (23 June 2022)), which is an application for automatic homology modeling [35], (iii) Phyre2 (http://www.sbg.bio.ic.ac.uk/phyre2/ (23 June 2022)), a remote homology detection service that creates structure as well as anticipates binding sites in the structure [36], (iv) MODWEB (https://modbase.compbio.ucsf.edu/modweb/ (23 June 2022)), a comparison modeling tool that utilizes PSI-BLAST to provide automated structural data [37], and (v) another similar platform, I-TASSER (https://zhanggroup.org/I-TASSER/ (23 June 2022)), an integrated platform for structure generation [38]. The most promising models from each tool were chosen for further evaluation.

2.6. Validation of ChrR Three-Dimensional Structure

Three reliable tools were employed to evaluate the models’ quality, including (i) ProQ (https://proq.bioinfo.se/ProQ/ (26 June 2022)), an algorithm based on neural networks for calculating the LGscore and MaxSub score of the structure [39], (ii) PROCHECK (http://www.ebi.ac.uk/thornton-srv/databases/pdbsum/Generate.html (26 June 2022)), which outlines the stereochemical characteristics of the structures [40], (iii) Verify Models (Dassault Systèmes BIOVIA, Discovery Studio Modeling Environment, Release 21.1, San Diego: Dassault Systèmes, 2021), a technique that calculates a DOPE score for each structure using the MODELER DOPE (Discrete Optimized Protein Energy) approach. A model with a lower score is statistically superior [41], and (iv) Verify Protein (Profiles-3D), which enables the evaluation of a protein sequence’s fitness in its current 3D environment [42]. The Verify score of an amino acid residue indicates whether or not the residue is located in the intended 3D environment. To obtain the most stable structure, the models performed energy minimization using the DS minimization algorithm. Using the molecular graphics technology of Discovery Studio Visualizer, the optimal structure produced after validation was subjected to a comprehensive structural analysis.

2.7. Molecular Docking and Interaction Studies

Considering the high sequence identity and structural similarities between the template (PDB ID: 1RTT_A) and the target protein, the active site was identified. To contribute in the localization of the binding site, MSA was used to explore the active site residues discussed in prior studies pertaining to ChrR [43,44]. Therefore, the predicted location was chosen as the most suitable binding site for ligand docking. To determine the optimal binding mechanism of the heavy metal ligand (Cr(VI)) to the receptor protein ChrR, docking studies were implemented by applying the CDOCKER tool from DS [45]. The preparation of the ligand and the receptor was completed by adding hydrogen, energy minimization, charge correction, and side chain refinement. According to studies, the metal-binding pocket in Gh-ChrR is close to the strongly bound FMN molecule [44]. The active site with x, y, and z coordinates 17.015425, −0.051868, and 66.484199, and a radius of 16.242105 Å were uploaded to define the active site. Using the View Interaction module from DS, close intramolecular contacts between Cr(VI) and ChrR that occurred between 2.5 and 3.5 Å were assessed to determine the stability of the docked complex. The binding energy of the receptor and the ligand was calculated using the calculate binding energy procedure:

Energy_Binding = Energy_Complex − Energy_Ligand − Energy_Receptor

2.8. In Silico Mutant Generation and Interaction with Cr(VI)

Eswaramoorth et al. described site-directed mutations in ChrR from E. coli that enhanced the chromate reductase activity [46]. Mutants were prepared by substituting active site residues to less bulky residue: Tyr128Asn, Glu146Thr, Tyr85Asn, and Arg185Met. Mutated models were generated and side chains were refined followed by energy minimization and stability assessment. The highest stable structure was docked with Cr(VI) for assessing binding affinity and intra-molecular interactions.

3. Results

3.1. Sequence Analysis

The amino acid sequence of bacterial chromate reductase from P. putida was retrieved from the UniProtKB database (https://www.uniprot.org/ (1 June 2022)) with the accession no. A0A1 × 0Z667. It consists of 186 amino acid residues. The UniProt analysis of ChrR from P. putida shows that it contains an FMN reductase domain from position 7–150.

Table 1. Physico-chemical parameters of chromate reductase.

Parameters	Measures
No. of amino acids	186
Molecular weight	20,282.20
Theoretical pI	8.53
Total number of negatively charged residues (Asp + Glu)	20
Total number of positively charged residues (Arg + Lys)	22
Ext. coefficient	25,565
Estimated half-life	30 h (mammalian reticulocytes, in vitro)
Instability index (II)	53.81
Aliphatic index	86.51
Grand average of hydropathicity (GRAVY)	−0.101

displays the calculated physicochemical parameters of the query sequences. The estimated instability index for A0A1X0Z667 is 40. ChrR’s virulence was predicted using the Virulentpred, MP3, and VICMpred prediction techniques. Consensus assessments for the sequence revealed it to be a non-pathogenic protein, as shown in

Table 2. Virulence factor prediction for ChrR from P. putida.

Tool	Prediction Approach	Prediction
VirulentPred	Amino Acid Composition-Based	0.4223 (Non-virulent)
	Dipeptide Composition-Based	−0.647 (Non-virulent)
	Higher-Order Dipeptide Composition-Based	−0.542
	Similarity-Based using PSI-BLAST	0 No hits obtained
	PSI-BLAST created PSSM Profiles	0.6117 (Virulent)
	Cascade of SVMs and PSI-BLAST	0.2162 (Non-virulent)
MP3	HMM	Non-pathogenic
	Hybrid	Non-pathogenic
	SVM	Pathogenic
VICMPred	Patterns + Compositions	−3.0806466 Metabolism molecule

.

3.2. Multiple Sequence Alignment and Evolutionary Analysis

In the current investigation, the eggNOG 5.0 program was implemented to identify the top seven orthologous sequences of the ChrR from P. putida. The orthologs were identified for ChrR, making them appropriate for phylogenetic analysis and MSA to identify conserved sites and reveal the evolutionary relationship, as shown in Figure 1. The Maximum Likelihood approach [47] and JTT matrix-based model [48] were utilized to deduce the ancestral states. The tree displays a range of probable amino acids (states) at each ancestral node based on their approximated likelihood at site 1. Only the most likely state for each node is displayed. The initial tree(s) for the heuristic search were constructed by applying the neighbor-joining and BioNJ algorithms to a matrix of pairwise distances determined using the JTT model, followed by identifying the topology with the highest log likelihood value. The rates were assumed to be the same across all sites (uniform rates option). The top seven orthologous amino acid sequences were selected for investigation, as shown in Figure 2.

3.3. Conserved Regions in the Sequence

The family and domain as assessed by CDD, SMART, Pfam, CDART, and InterProScan are listed in

Table 3. Identified domains in ChrR.

Tool	Domain	Residues
SMART	FMN_red	5–153
	Flavodxin_2	6–143
Pfam	FMN_red	5–153
CDD	SsuE (NAD(P)H-dependent FMN reductase)	6–177
CDART	SsuE (NAD(P)H-dependent FMN reductase)	3–188
InterProScan	FMN_rdtase-like	7–150
	FMN_red	7–150

. ChrR was found to possess a flavine mononucleotide (FMN) reductase domain in the region 5–153 and a flavodoxin domain in the region 6–143. The identified motifs and their locations using MEME-Suite are: PEYNRDV (17–23), PEYNRDV (46–52), WGGFGEEI (61–68), PEYNRDV (79–85), WGGFGEEI (120–127), KDRPCLQ (136–142), and KDRPCLQ (164–170).

3.4. Secondary Structure Content in the Sequence

As illustrated in Figure 3, the anticipated secondary structure of ChrR from P. putida consists of 41.39% helices, 16.12% sheets, and 42.48% loops. This is supported by the published 3D structure (PDB ID: 3U7R) and the topology of the sequence with the highest degree of similarity.

3.5. Structure of ChrR from P. putida

Using PSI-BLAST to look for similar sequences in the PDB database, it was found that the A0A1X0Z667 protein sequence is most similar to ChrR from Paracoccus denitrificans (PDB ID:3U7R_A), an NADPH-dependent FMN reductase with the highest similarity of 48.08% and maximum query coverage of 95%. Another sequence with second highest similarity of 43% and sequence coverage of 97%, which is a putative NADH-dependent reductase belonging to Pseudomonas aeruginosa, was considered as a template sequence for the homology modeling of ChrR P. putida. Five structures of ChrR from P. putida were produced using Phyre2, ModWeb, SwissModel, Discovery Studio, and I-Tasser. The optimal model created by each tool was validated using ProQ, Verify Protein (MODELER and Profile3D), and PROCHECK. The model generated by I-Tasser contained all 186 residues and the lowest MODELER score; however, it showed only 75.5% of residues in the most favored region of the Ramachandran Plot. The assessment thus carried out in every aspect for each of the generated model revealed the Phyre2 model as the best model for further investigations, as presented in

Table 4. Validation of homology-modeled structures of ChrR from P. putida.

Modeling Tool	Residues	ProQ		Verify Protein		PROCHECK
		LG	Max Sub	MODELER	Profile 3D	Most Favored	Additionally Allowed	Generous Allowed	Dis Allowed Region
Phyre2	4–183	−0.31	8.58	−20,023.46	87.74	87.6	9.8	2.6	0.0
ModWeb	4–183	−0.31	8.2	−20,425.416	87.97	91.5	7.8	0.7	0.0
SwissModel	1–186	−0.36	7.67	−20,981.69	86.28	92.5	4.4	0.6	1.3
Discovery Studio	185	−0.604	8.5	−19,997.08	80.85	77.25	20.3	1.3	1.3
I-Tasser	1–186	−0.55	8.28	−21,941.75	84.23	75.5	22.0	2.5	0.0

.

The homology-modeled structure comprised 180 residues that start from Val4 and end at Asn183. The model contained six α-helices and five β-sheets connected by seven loops. The structure was found to have a flavodoxin fold. Moreover, CATH characterized it to belong to the Alpha Beta class, 3-Layer(aba) Sandwich, Rossmann fold, and Flavodoxin domain. In addition, according to InterProScan, it belongs to the Flavoprotein-like sf superfamily. It was observed that all four β-sheets are arranged parallelly at the center of the structure and surrounded by α2, α3, and α4 helices at one side, and on the other side lies α1 and α5. The cavity formed near the c-terminus end of the β-sheet by Loop1, Loop2, and Loop3, corresponds to the FMN-binding site. A nucleotide-binding motif “GSLRKESYN” is located at Gly12 to Asn20 on Loop1. ChrR belongs to the putative NADH_dh2 family of flavin-binding quinone reductases, and it carries the signature sequence motif “LFVTPEYNXXXXXXXXLKNAIDXXS” at amino acid positions 75–97 [49], as shown in Figure 4. Cr(VI) is bound in the binding pocket of ChrR with a -CDocker score of 27.0658, -CDocker Interaction Energy of 11.6155, and binding energy of 0.00 kcal/mol. Cr(VI) was found to interact with the O atom of Glu80 at a distance of 1.55 Å, as illustrated in Figure 5.

3.6. In Silico Site-Directed Mutagenesis

Corresponding residues to the reported mutant residues: Tyr85, Arg125, Tyr128, and Glu146, were identified through sequence alignment [46]. The corresponding residues are Arg83, Gly124, His127, and Glu145, respectively. These residues were mutated to the highest stable mutant. To assess the interaction, the mutant models were docked with Cr(VI). Arg83, when mutated to Trp, showed a stabilizing effect with a mutation energy of −1.14 kcal/mol, as shown in Table S1. Gly124, when mutated to Ile, Trp, Phe, Val, Leu, Tyr, Gln, Glu, His, Asp, Thr, and Ala, showed a stabilizing effect. However, the lowest mutation energy was observed in Gly124Ile at −2.77 kcal/mol, as shown in Table S2. The stabilizing effect of the His127 mutation to Trp, Phe, Ile, Leu, Cys, Val, Arg, Lys, Gln, Glu, Tyr, Asn, Asp, and Lys is shown in Table S3. However, the lowest mutation energy, −1.62 kcal/mol, is observed in mutant His127Trp. The mutation of Glu145 to other amino acids showed a neutral and destabilizing effect. Therefore, the other stabilizing mutants were analyzed for interaction with Cr(VI), as shown in Table S4. Mutant models of Arg83Trp, Gly124Ile, and His127Trp were generated using the “Build Mutant” protocol from DS. All of the mutant structures were energy minimized. The one with the lowest energy was chosen among the possible docking models for Cr(VI). As seen in Figure 6, all of the mutants were shown to participate in close intra-molecular contact with the O atom of Glu80 at a distance of only 1.15 Å.

4. Discussion

The sequence analysis revealed an instability index greater than 40, implying it to be an unstable protein. Consensus assessments for the pathogenicity of the enzyme sequence confirmed that it was non-pathogenic and safe for human usage. It was noted from the MSA that amino acids belonging to the signature motif are mostly conserved among the orthologs. The phylogenetic tree derived from this research demonstrated that ChrR from Paracoccus denitrificans was closely linked to P. putida. It is generally believed that amongst orthologs, there is a high degree of conservation in the functional specialization of proteins. Based on this hypothesis, we were able to locate the amino acids that determine the specificity of protein–ligand interactions. Identifying these residues is critical for studying molecular recognition mechanisms and designing effective proteins and drugs [47]. According to the phylogenetic tree, ChrR was discovered to be closely associated with NADPH-dependent FMN reductase from the bacteria P. denitrificans. In addition, FMN-dependent NADPH-quinone reductase (azoR) from B. cohnii was found to be the most probable ancestor. ChrR was found to possess a flavine mononucleotide (FMN) reductase domain in regions 5–153 and a flavodoxin domain in regions 6–143. The FMN reductase reduces flavins by utilizing NADH or NADPH. It is engaged in several biological activities, such as bacterial bioluminescence, that involve facilitated flavin for specific functions. The flavodoxin is composed of structural domains with a flavodoxin-like fold, which is present in various proteins, including flavodoxins, and FMN-dependent NADH-azo compound oxidoreductases. The flavodoxin fold is defined by an open twisted/alpha beta structure composed of five parallel beta-sheets connected by encircling alpha helices. Flavodoxins are proteins that transmit electrons and work in various electron transport systems. They are interchangeable with ferredoxins in their ability to bind one FMN molecule, which acts as a redox-active prosthetic group [50]. Therefore, the presence of an FMN-binding domain assures binding and reduction of Cr in the presence of FMN. The sequence contains Rossman fold topology; it is one of the most prevalent and extensively found super-secondary structures. The hydrogen-bonded beta strands form a beta-sheet, and the alpha-helical segments alternate with the beta strand segments to create the structure. The portion of the Rossman fold most preserved is the initial beta–alpha–beta fold. This region is also known as an “ADP-binding fold” as it touches the ADP portion of dinucleotides, including FAD, NAD, and NADP [51]. A thorough investigation of the ChrR sequence revealed its function in metal ion binding, iron–sulfur cluster binding, oxidoreductase activity, and catalytic activity.

The secondary structure prediction is generally accurate and considerably simpler to solve than the three-dimensional structure prediction. It is essential for the precise prediction of 3D structures. Furthermore, secondary structure prediction can potentially evaluate the reliability of a model constructed using a (tertiary) structure prediction technique [52]. The secondary structure content of the modeled protein revealed a higher range of helices than sheets, indicating considerable compactness in the structure.

The docking studies highlight the substrate specificity of flavoenzymes, which is often broad. They have been predominantly described in the context of alternative substrates such as nitroreductase NfsA, ferric reductase FerB, and numerous chromate reductases [49,53]. Ackerley et al. reported a quantitative transformation of chromate to Cr(III) by ChrR, in which the final product of ChrR-catalyzed reactions was Cr(III); it is evident that the enzyme catalyzed more than a simple one-electron reduction of Cr(VI) to Cr(V) [54]. Therefore, the molecular docking approach assessed binding and intra-molecular interactions of Cr(VI). Most electron donors for reducing Cr(VI) are organic molecules. However, inorganic compounds have also been found to be reducing equivalents. Significant electron donors for Cr(VI) reduction include fructose, lactose, lactate, glucose, pyruvate, citrate, acetate, glycerol, formate, reduced glutathione, NADH/NADPH, and others. Gluconacetobacter hansenii (Gh)-ChrR is an FMN-containing chromate reductase that undergoes structural rearrangement in response to the interaction of chromate anion and NADH. With this rearrangement, both species may attach simultaneously for practical enzyme cycling; otherwise, the electron donor would have occupied the binding site. In addition, Cr(VI) reduction by microorganisms using glucose as the donor material is quite popular. The present study reveals the binding of Cr (VI) in the binding pocket and forming metal interactions with Glu80, a crucial residue in the signature motif for FMN binding.

Engineered proteins, particularly enzymes, are being used more frequently in various industries due to their catalytic ability, selective ligand binding, and capabilities as materials and food additives. The urge to engineer or generate proteins with higher stability, activity, and specificity has increased along with the number of possible applications for engineered proteins. As the application of protein technology develops, exploiting the potential advantages of modulating remote regions will become imperative [55]. Mutations with enhanced chromate reductase activity relative to the wild-type reported by Eswaramoorthy et al. were compared with ChrR from P. putida.

Bioremediation strategies, such as biostimulation, is a potential bioremediation technology [56]. In this method, nutrients are given to the surrounding environment, including aquifers, to promote the growth of microorganisms. The change improves bioremediation, but excessive biomass may clog subsurface pores, reducing the effectiveness of the cleanup. Mixed garbage hinders remediating bacteria and enzymes. Hence, the biostimulation of such areas is ineffective. Toxic intermediates formed during chromate reduction are harmful to remediate microorganisms, complicating chromate clean-up. One potential remedy for these problems is genetic and protein engineering techniques [57]. Slow-growing bacteria can have their desired genes expressed to their full potential with particular promoters to reduce biomass accumulation and clog. Enhanced enzymes that reduce chromate more effectively, with less toxicity to the bacteria that are doing the remediation, and the ability to function in the presence of additional contaminants can also be produced by the protein engineering of bacterial chromate reductases [58].

5. Conclusions

The complicated issue of reducing Cr(VI) contamination from the environment has contributed to development of a range of bioremediation techniques, particularly microorganism-based reduction techniques. Pseudomonas putida is an adaptive bacterium. Due to its robust metabolism, tolerance to toxic elements and oxidative stress, and flexibility of genetic modification, this bacterium is becoming a cell factory, manufacturing natural products with varied biological purposes. The present study provides a validated three-dimensional structure of ChrR from P. putida and mutant models, stabilizing the system for improved binding with Cr(VI) to reduce it to Cr(III). Indeed, the evidence reported here points to a crucial role of Glu80 in binding with Cr(VI). We have shown that the mutant models Arg83Trp, Gly124Ile, and His127Trp retain their property of binding with Cr(VI) for reducing it to Cr(III). The wild-type ChrR from P. putida was observed to be unstable. However, the mutants stabilized its structure and retained the metallic interaction formation between the critical amino acid residue of the signature motif. These mutants can be further validated for protein engineering for introducing stable ChrR that can survive and reduce Cr(VI) in the harsh climatic conditions of Madinah, Saudi Arabia. Our observations revealed that ChrR from P. putida is a potential enzyme for the bioengineering and bioremediation of Cr; our future research will focus on the role of other non-pathogenic microorganisms in the bioremediation of heavy metals contaminating the groundwater of Madinah.