A lignite sample from the coal mine in Neyveli, India (11°35′34.44″N and 79°29′29.04″E), was collected. The metagenomic DNA was isolated from the lignite sample using the PowerMax soil DNeasy kit (QIAGEN). The isolated metagenomic DNA was then subjected to shotgun sequencing on an Illumina HiSeqTM 2000 platform to generate paired-end sequences.

The forward and reverse end reads in the FASTQ format were used as the input in the FASTQC tool (Andrews, 2010). The generated output HTML files were merged using MultiQC (Ewels et al., 2016) to create a single HTML file report containing the quality statistics of the reads. The forward and reverse FASTQ files are the input for the MEGAHIT assembler (Li et al., 2015). A k-mer value of 99 and a minimum contig length of 200 parameters were assigned. The output was generated as contigs in a single FASTA file. The obtained contigs were analysed in the QUAST tool (Gurevich et al., 2013) for the number and size of the contigs.

The contigs in FASTA format were annotated using the SqueezeMeta tool, an automated pipeline (Tamames and Puente-Sánchez, 2019). First, protein-coding genes were predicted from the contigs using Prodigal, and amino acid and nucleotide sequences were generated in the FASTA files. The results of these annotated nucleotide sequences were automatically loaded as input into Diamond, which searched the GenBank nr database for taxonomical classification and the KEGG database for functional annotation. The term “Nitrilase” was searched using the grep script in tab-separated value files. KEGG IDs and contig IDs were noted for identifying taxonomy and extracting the nitrilase nucleotide and amino acid sequences.

The identified amino acid sequence in FASTA format was uploaded in BlastP (Johnson et al., 2008) and ran against the NCBI database of protein reference sequences to find similar sequences. Similar sequences were chosen based on the >50% identity of the matches with the query sequence. These sequences above the threshold were downloaded in FASTA format in a single file along with the query sequence. This FASTA file was uploaded using MEGA 11 software (Tamura et al., 2021) and was aligned using Clustal Omega. The evolutionary history of these sequences was created using the neighbour-joining method by selecting the phylogeny tab in the software application.

The amino acid sequence in FASTA format was submitted with the default parameters in the NCBI–CDD (Marchler-Bauer et al., 2017) and InterProScan (Quevillon et al., 2005) databases to predict the homologous superfamily, conserved domain, conserved region, Gene Ontology, and NCBI-CDD from the query sequence of amino acids.

The amino acid sequence was pasted in FASTA format into a query box of the ExPASy’s ProtParam (Gasteiger et al., 2003) server and submitted to identify the physical and chemical properties of the protein sequence. The server page quantifies the number of amino acids, the molecular weight, the number of negatively charged residues, the instability index, the theoretical pI, and the grand average of hydropathicity.

The secondary structure was predicted by uploading an amino acid sequence in FASTA format to the NetSurfP tool with the default parameters (Høie et al., 2022). The amino acid sequences were used in the software package Chimera X version 1.4 (Pettersen et al., 2021; Goddard et al., 2018) to perform the AlphaFold2 tool for 3D structure prediction. This was predicted in the software programme by choosing the AlphaFold2 option from the Structure Prediction tab under Tools. The amino acid sequence was uploaded in the query box and submitted with the default parameters. The predicted tertiary structure was uploaded into the PROCHECK tool (Laskowski et al., 1993) to create the Ramachandran plot, which checks to validate the stereochemical quality of the protein structure.

The predicted structure was submitted for molecular dynamics simulation on the WebGRO server (Bekker et al., 1993) to check its stability. Using the GROMOS96 43a1 force field settings, the complex system was solvated using a simple point charge (SPC) water model in a triclinic periodic box. The complex system was maintained at a salt concentration of 0.15 M by adding a suitable amount of Na⁺ and Cl⁻ counterions. Using the steepest descent approach, energy reduction was achieved in 5,000 steps. Constant amount, volume temperature (NVT/NPT), and pressure equilibrium types were used. The temperature was set to 300 K, and the pressure was set to 1.0 bar. The simulation time was 100 ns and was conducted with 1,000 frames per simulation. Finally, the simulation result was analysed based on the time-dependent root mean square deviation (RMSD) of the given structure and the root mean square fluctuation (RMSF) of each residue.

The homepage of the Protein Data Bank (PDB) (Burley et al., 2017) was searched for the 3D nitrilase protein structure. X-ray diffraction was selected using the filter option, and the crystal structure of nitrilase proteins was retrieved in PDB format. Nitriles were downloaded from the PubChem database (Kim et al., 2020) in SDF format. The ligands were converted to the PDB format using PyMOL (Schrödinger, 2000), and protein and ligand formats were changed to the PDBQT format for molecular docking using the AutoDockTools-1.5.7 tool.

The protein structures’ active site residues were predicted by uploading the PDB files of nitrilase proteins in the computed atlas of the protein surface topography −3.0 (CASTp) server (Tian et al., 2018). The alpha shape theory’s pocket algorithm calculated the active pockets or binding sites, with large pockets with high volumes likely to contain enzyme-binding sites for the interaction between proteins and ligands.

In AutoDockTools-1.5.7, the protein structures in the PDB format were used as input to create a grid file by placing a grid box in the protein’s predicted active site for the binding of ligand molecules. The inputs of a protein, a ligand in PDBQT format, and a grid file in text format were used to perform a molecular docking analysis in AutoDock Vina (Eberhardt et al., 2021) with the default parameters. The output was the protein–ligand complex in PDB format. The complex in PDB format was used to study the interaction of amino acids with ligands using the LigPlot + tool (Laskowski and Swindells, 2011).

FASTQC and MULTIQC determined the paired-end reads to be within the Phred score value of 36, considering that the raw reads to be of standard quality (Bin Kwong et al., 2017). The quality reads of 150-bp length from 32 GB of data were processed for assembly. The statistics of the assembled contigs were evaluated using the QUAST tool, resulting in 226 contigs with >50,000 bp, 1,240 contigs with >25,000 bp, 5,265 contigs with >10,000 bp, 11,596 contigs with >5000 bp, and 73,144 contigs with >1000 bp.

The phyla observed in significant numbers belong to Proteobacteria (76%), Actinobacteria (8%), Firmicutes (8%), Spirochaetes (2%), Bacteroidetes (1.5%), Chloroflexi (0.8%), Planctomycetes (0.5%), Cyanobacteria (0.2%), [Thermi] (0.1%), Fusobacteria (0.1%), and Acidobacteria (0.1%). Using “grep,” the results from SqueezeMeta were searched for nitrilase, and then organisms involved in nitrilase enzymes were identified. In total, 27 organisms were identified, of which four were from unclassified microorganisms and belonged to the KEGG ID K01501 with the metabolic pathway number EC:3.5.5.1. For the novel protein approach, long amino acid sequences and unclassified microorganisms were preferred. A single amino acid sequence of 261 base pairs encoding an unclassified Alphaproteobacteria was selected.

Using BlastP, the target protein sequence was compared to the NCBI database of protein reference sequences. Eight sequences above the threshold were downloaded from NCBI-BLAST, and a phylogenetic tree was built along with the identified nitrilase sequence using the MEGA11 tool. Figure 2 shows the multiple sequence alignment of all the sequences. The branch lengths marked next to the branches on the tree are shown to scale and are in the same units as the evolutionary distances, which are used to estimate the phylogenetic tree (Figure 3). The scale value of 0.050 shows the genetic change between the protein sequences. The branch of the nitrilase coal metagenome is the longest, with a length of 0.30, and is therefore known to have greater genetic change than other protein sequences. The software application calculated evolutionary distances using the Poisson correction method and amino acid substitutions per site. This analysis involved nine protein sequences. There were a total of 309 positions in the final data set, and 139 residues were conserved among the protein sequences.

According to the analysis of the InterProScan database and NCBI-CDD results, a carbon–nitrogen hydrolase domain was found in the 4–19 amino acid position in the protein sequence, indicating it to be a member of the nitrilase superfamily. InterProScan classifies three domains: cellular component, molecular function, and biological process. The predicted biological process for the identified protein is involved in the metabolic process of the nitrogen compounds (GO: 0006807), and the molecular function is involved in catalytic activity (GO:0003824). Based on the unclassification, amino acid sequence length, conserved region, and domain analysis, it is considered to be a novel nitrilase enzyme.

Using the ProtParam tool, the physiochemical characteristics of the novel protein sequence were calculated. The compositions of the amino acids are shown in Figure 4. The result shows that the negatively charged aspartic acid and glutamic acid were higher than the positively charged arginine and lysine. The protein has an acidic nature because of its 4.86 isoelectric points. The total number of atoms is around 4,010, including 1,275 carbon, 1,990 hydrogen, 354 nitrogen, 385 oxygen, and 6 sulphur atoms. Proteins with a GRAVY score below zero are hydrophilic, while proteins with a GRAVY score above zero are hydrophobic. The novel protein has a GRAVY value of −0.126, which indicates that the protein is inherently hydrophilic. The protein has an instability index of 53.02, which indicates that it is unstable because a protein with an instability index of 40 or less was stable in the test tube. The sequence of amino acids was also used to estimate the protein’s half-life. For yeast cells, the half-life period is > 20 h (in vivo), and for E. coli, the half-life period is >10 h (in vivo) (Gasteiger et al., 2005).

The results from the NetSurfP version 2.0 server revealed that the nitrilase protein has 10 helices, 16 strands, and 23 coils in its secondary structure (Figure 5).

Using the amino acid sequence with the default parameters, the nitrilase protein structure was predicted in AlphaFold2 (Figure 6). In AlphaFold2, the predicted local distance difference test (pLDDT) calculates the mean confidence value ranging from 0 to 100. Here, the predicted nitrilase’s mean pLDDT score was 95.8. A mean confidence value of 90 and above is said to agree with an experimental structure (Jumper et al., 2021). The Ramachandran plot was utilised in the PROCHECK tool, and the modelled protein was evaluated and validated (Supplementary Material S1). In the model, 92.7% of residues were found in the most favoured regions (A, B, and L), 6.8% in the additional allowed regions (a, b, l, and p), 0.5% in the generously allowed regions (∼a, ∼b, ∼l, and ∼p), and no residues in the disallowed regions. As a result, the percentage distribution of amino acid residues revealed that the predicted nitrilase structure is of high quality.

To analyse the flexibility and stability of the best-predicted protein structure provided by AlphaFold2, a time-dependent molecular dynamic simulation was conducted at 100 ns employing the GROMACS forcefield on the WebGRO server. The RMSD value of the nitrilase protein structure showed equilibrium with an average of 0.2 nm. The largest oscillation in RMSD was observed in the 0–10 ns range. Afterwards, the 10-ns RMSD value was stabilised up to 100 ns, with an average value of 0.5 nm (Figure 7). The RMSF value was then assessed to analyse the structural flexibility of the atoms in the backbone of the proteins. The obtained data showed that there were fluctuations in the residues present in the loops of the protein structure (RMSF ≤ 0.5 nm), and this indicates that the complex was flexible in these loop regions (Supplementary Material S2). The conformational changes in the loop structure caused by the flexibility have no impact on the protein structure (Li et al., 2011).

The predicted nitrilase protein and the retrieved crystal structures of 3WUY (Zhang et al., 2014), 3IW3, and 3IVZ (Raczynska et al., 2011) were used to check the binding scores individually by molecular docking analysis and compared. The three structures were determined by X-ray diffraction analysis. The protein structure (3WUY), which was identified from Synechocystis sp. PCC 6803 substr. Kazusa, comprises 349 amino acids and has a resolution of 3.10 Å. The other two protein structures, 3IW3 and 3IVZ, were identified in Pyrococcus abyssi GE5 and have a length of 262 amino acid sequences with resolutions of 1.80 Å and 1.57 Å, respectively. For the protein–ligand interaction study, nitriles such as acrylonitrile, benzonitrile, dichlobenil, fumaronitrile, malanonitrile, and succinonitrile were retrieved. Proteins and ligand molecules were converted to the PDBQT format and were ready for docking analysis.

The CASTp server identified the active site for the predicted nitrilase protein and the crystal structure of the proteins. The surface area measurement and cavity volumes were predicted. For the predicted protein, the area of the active site was 677.110 Å2, and the volume was 645.046 Å3. The area of the active site for 3WUY, 3IW3, and 3IVZ was 983.426 Å2, 61.447 Å2, and 188.905 Å2 and the volume was 1,129.287 Å3, 31.456 Å3, and 60.553 Å3, respectively. All the pockets chosen for the proteins have a high surface area and volume for the specific enzyme-binding site (Barman et al., 2020).

Using AutoDockTools-1.5.7, a grid file was generated for all the proteins. When the model overlapped with the template, an acceptable range of RMSD was determined to be 2.0, which is regarded as satisfactory docking (Castro-Alvarez et al., 2017). The docking analysis was completed using AutoDock Vina, as tabulated in Table 1. The docking results produced nine ligand-binding poses, of which the one with the lowest RMSD (value 0) was selected, indicating a true binding pose. The protein–ligand interaction showed that these proteins’ docking scores are almost similar when they bind to six different ligands.

The protein–ligand complex of the predicted nitrilase protein was analysed to determine the nature of the interactions between the amino acids and the ligands using LigPlot + (Figure 8).

Nitrile (C≡N) conversion involves a nucleophilic attack on the electrophilic carbon (electrophilic) by the nucleophilic group present in the side chain of amino acids in the active site, resulting in the formation of carboxyl groups (COOH). Cysteine, serine, threonine, tyrosine, glutamic acid, aspartic acid, lysine, arginine, and histidine have nucleophilic R groups and can act as nucleophilic donors. Arginine and tyrosine amino acids serve as functional converters in benzonitrile and succinonitrile. Threonine and glutamic acid play a functional role in the catalysis of fumaronitrile and serine in malononitrile. The nucleophilic attack on acrylonitrile and dichlobenil is mediated by histidine and threonine, respectively.