A study on transcriptome analysis of G. max by, Wang et al., 2014 revealed the induction of hundreds of genes and transcript diversity in response to S. litura feeding. Out of diverse induced transcripts, we have selected biotic stress-responsive differentially expressed CYP450 genes to analyze their role in G. max- S. litura interaction. Protein and nucleotide sequences of selected 16 GmCYPs were fetched from the Phytozome database version 13 (https://phytozome-next.jgi.doe.gov/) and also manually evaluated using NCBI-BLAST tool (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Molecular weight (MW), isoelectric point (pI), instability index, and hydrophilicity (GRAVY value) of candidate sequences were computed employing the ProtParam tool (https://web.expasy.org/protparam/) of ExPASy (https://www.expasy.org/). ExPASy server UniProtKB/Swiss-Prot (https://www.expasy.org/resources/uniprotkb-swiss-prot) was used to calculate the length of amino acids for all the listed GmCYPs sequences (Gasteiger et al., 2005). ExPasy-Prosite server was used to identify functional sites. Deep-Loc Tool (https://services.healthtech.dtu.dk/service.php?DeepLoc-1.0) was employed to analyze the subcellular localization of listed GmCYPs to gain insight into their functional mechanisms. Subcellular localization of all GmCYPs was represented as a heatmap using TBtool software. Subsequently, sequences were uploaded on the TMHMM version 2.0 online tool (https://services.healthtech.dtu.dk/service.php?TMHMM-2.0) for trans-membrane domain analysis (Wang et al., 2022).

To evaluate the evolutionary connections among leguminous dicots, model plants, monocots, homologous protein sequences of listed GmCYPs for Arabidopsis thaliana, Nicotiana attenuata, Medicago truncatula and Oryza sativa were acquired from NCBI database using BLAST search (https://www.ncbi.nlm.nih.gov/guide/proteins/). Subsequently, the Clustal omega tool was used to align amino acid sequences (https://www.ebi.ac.uk/Tools/msa/clustalo/). The maximum-likelihood phylogenetic tree was produced based on the result of multiple sequence alignments using MEGA 7.0 (Molecular Evolutionary Genetic Analysis) software (https://megasoftware.net/) with bootstrap replication value 1000, P-distance mode and partial deletion (Kumar et al., 2016). Finally, the phylogenetic tree was embellished, implying iTOL (https://itol.embl.de/upload.cgi) tool (Letunic and Bork, 2007).

MCScanX software was employed to identify duplicated CYP gene pairs found in the genomes of G. max, A. thaliana, M. truncatula, N. atteneuata and O. sativa through their non-redundant pairing, thereby generating five pairs (G. max/G. max, G. max/A. thaliana, G. max/M. truncatula, G. max/N. atteneuata and G. max/O. sativa) (Wang et al., 2012). The software employed an alignment file (generated by executing the all-versus-all BLASTP program for each of five pairs with e-value 10-5) and genome annotation data (retrieved from RefSeq of NCBI database) as input to generate a collinearity file. Collinearity file, which includes collinear pairs found in CYP genes for each species pair, was utilised to map duplicate gene pairs through the Circos software (http://circos.ca/software/) package (Krzywinski et al., 2009). To investigate selection pressure employed in GmCYPs duplication, the Ka/Ks ratio (synonymous/non-synonymous substitution ratio) was calculated through PAL2NAL software (http://www.bork.embl.de/pal2nal/) (Suyama et al., 2006).

Exon-intron organization of listed GmCYPs was illustrated using GSDS (Gene Structure Display Server) (http://gsds.gao-lab.org/), where genomic sequences of all the genes were aligned with their corresponding cDNA sequences. TBtool software was used to visualize the gene structure of 16 GmCYPs (Hu et al., 2015). Listed GmCYPs were mapped to soybean chromosomes based on the generic feature format version file (GFF) retrieved from the Phytozome database. Chromosomal position and relative distances of genes were spotted using Ritchie lab-Phenogram software (http://visualization.ritchielab.org/phenograms/plot) (Abbas et al., 2022).

Functional domain in GmCYPs was searched by submitting protein sequences to Pfam: Protein families database tool (http://pfam-legacy.xfam.org/). InterPro 91.0 database (https://www.ebi.ac.uk/interpro/) and NCBI Conserved domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) were also used for the detection of the conserved domains (Mistry et al., 2021). Later, MEME: Multiple Em for Motif Elicitation software (https://meme-suite.org/meme/doc/meme.html) was employed to predict conserved motifs and generate logo plots of listed GmCYPs sequences selecting the number of motifs: 10 and rest parameters as default (Bailey et al., 2009). Multiple sequence alignment (MSA) with amino acid sequences of candidate GmCYPs was performed to investigate existing conserved motifs utilizing Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/).

The 1000 bp nucleotide sequences upstream of translational start codon of each listed GmCYPs were retrieved from Phytozome v13 database. The cis-acting elements in listed GmCYPs sequences were predicted employing the Plant Care database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). Further, the results of cis-acting element was analysed using TBtool software (Zhang et al., 2022). psRNA Target tool (https://www.zhaolab.org/psRNATarget/; Lescot et al., 2002) was employed for anticipation of miRNA that may target and regulate GmCYPs expression with default parameter settings (Dai and Zhao, 2011). The mature sequences of miRNA of G. max were accessed from the miRBase i.e. miRNA database (https://mirbase.org/). The GmCYPs mRNA sequences and miRNA sequences were submitted to psRNA Target analysis server to characterize GmCYP’s miRNA target sites (Chen et al., 2020, Dai and Zhao, 2011).

Transcription factor (TF) binding sites prediction in the GmCYPs promoter region was executed using PlantTFDB online server (http://planttfdb.gao-lab.org/); using default parameters and p-value of ≤ 1e-6 (Malviya et al., 2023).

To analyse the organ-specific expression profiles of genes, RNAseq data related to differential expression in different tissues was retrieved from SoyBase Expression Explorer for each GmCYPs (https://www.soybase.org/expression). Expression values were measured as Reads/Kb/Million (RPKM)) for various organs, including the seed, stem, nodule, flower, young leaf, and green pod (Waese et al., 2017). Normalized RPKM values of the listed GmCYPs were used to generate a heat map using TBtool.

The protein-protein interaction network of GmCYPs was deduced using the String database (https://string-db.org/). For PPI, G. max was selected as the model plant, with medium confidence = 0.400 and the remaining parameters at the default settings (Szklarczyk et al., 2021).

Foliar application of an equal volume of 50 mM ethephon, 5 mM salicylic acid, 100 mM of methyl-jasmonate, and 0.5% (v/v) ethanol in water/only water was administrated on one-month-old healthy soybean plants kept in individual enclosures under identical conditions (Singh et al., 2008; Yadav et al., 2022). The treated and control samples were collected after 24 hours and immediately frozen at -80°C. The samples were collected in biological triplicates for subsequent analysis. Data for gene expression analysis was interpreted employing the 2^-ΔΔCT method. T-test (P < 0.05) and ANOVA were applied for statistical analysis, while differentially expressed genes were selected based on a criterion of more than two-fold, either induction or reduction.

Three-dimensional structures of GmCYPs were predicted using Phyre2, Swiss modeller and i- TASSER by submitting protein sequences. Structures were analyzed based on different parameters such as query coverage, template size and origin, stability of the predicted proteins and validated by Ramachandran plot (Kelley et al., 2015; Waterhouse et al., 2018; Zhang, 2008). After which, the best model were utilized for further analysis such as identification of conserved motifs and highlighting them on the most accurately modelled structures using Pymol (DeLano, 2002). Docking of heme prosthetic factor with one of the randomly selected GmCYP was executed using Maestro Glide docking software (Schrodinger’s tool) (Bhachoo and Beuming, 2017). Prior to docking, ligand (Heme) and GmCYP structures were optimized by Maestro’s Ligprep and Protein Preparation wizard, respectively. Receptor grid generation was constructed by selecting residues of heme-binding motifs. Arabidopsis CYP94C1 was taken as a reference for this analysis.

Information about substrates of GmCYPs was retrieved from literature and databases survey, and 3-D structures of GmCYPs substrates were fetched from PubChem. To show the interaction between GmCYPs and their substrates, the InChl format of the substrates fetched from PubChem was converted to. pdb format employing Open Babel-Chemical file format (http://www.cheminfo.org/Chemistry/Cheminformatics/FormatConverter/index.html) (O’Boyle et al., 2011; Kim et al., 2019). Patch Dock was applied to perform interaction studies with parameters: RMSD value= 4 and protein-small ligand complex type, and 100 top results. The ten highest-rank Patch dock docking outcomes were refined through FireDock (Mashiach et al., 2008; Schneidman-Duhovny et al., 2005).

Gene IDs for all S. litura-induced GmCYPs were accessible from NCBI and phytozome database, from where their protein and nucleotide sequences were retrieved. A self-BLAST search was also conducted to eliminate redundant sequences that finally identified 16 of the S. litura-inducible GmCYP450s involved in biotic stress response. The naming of genes was based on their homology to Arabidopsis CYP450s and Phytozome gene description. Detailed information of 16 candidate GmCYPs, including gene ID, Phytozome ID, Uniprot ID, and the number of peptide sequences along with physiochemical properties, are represented in Table 1 . The candidate GmCYPs range from 6152 bp (GmCYP78A5) to 2164 bp (GmCYP78A3), and CDS sequences ranged between 1398 bp (GmCYP85A1) to 1746 bp (GmCYP78A5). The full length of peptide sequences varied from 466 (GmCYP85A1) to 582 (GmCYP78A5) amino peptides with molecular weight (MW) varying from 53.375 (GmCYP85A1) to 65.343 (GmCYP78A5) kDa and their isoelectric point (pI) values varying from 6.54 (GmCYP82A2) to 9.27 (GmCYP90A1), indicating that maximum proteins were basic. However, instability index values indicate that most GmCYPs are stable (Instability index > 40) except GmCYP71AU50, GmCYP93A3, GmCYP89A2, and GmCYP71D10. In addition, 13 genes were hydrophilic (GRAVY value less than zero), while GmCYP78A5-like, GmCYP76C1-like Determination of subcellular localization were predicted using reticulum (ER). Their presence was also predicted in the chloroplast, mitochondria, cell membrane, vacuole, golgi apparatus and nucleus with minimal expression, as depicted in Supplementary Figure 1 .

A phylogenetic tree was constructed from G. max CYP450s genes (16) and their homologous members from A. thaliana (14), M. truncatula (14), O. sativa (14) and N. tabacum (15); to understand the evolutionary relationship among them ( Figure 1A ). The findings show that CYPs have evolved dramatically throughout the time series, as seen by distinct clades in tree, implying that CYPs play a variety of functions in various plant systems. These proteins were categorized into 12 clades; among them, Clan 71 has a maximum of 45 members, and Clan 80, Clan 72, Clan 92, and CYP736 contained 10, 4, 11 and 3 members, respectively. GmCYPs were dispersed across the phylogenetic tree and belonged to different clades illustrating evolutionary relatedness with CYPs from A. thaliana, M. truncatula, O. sativa and N. tabacum. However, the maximum like-hood tree indicated that G. max and M. truncatula are closely related evolutionarily (having the highest bootstrap values), inferring they are members of the same family. Gene duplications are crucial for expanding and evolving gene families, driving species divergence and new gene functions (Innan and Kondrashov, 2010). Examination of gene duplication occurrence in the G. max genome revealed that 15 pairs of GmCYPs underwent segmental duplication i.e. duplicated genes were found on different chromosomes. A total of 6 pairs of duplicated genes between GmCYPs and genes from A. thaliana while 12 pairs of duplicated genes between GmCYPs and genes from M. truncatula were observed ( Figure 1B ). The selection pressure was interpreted by calculating Ka/Ks which varied from 0.04918 (NP_850337.1/XP_003521103.1) to 0.332959 (XP_006582383.1/XP_039688169.1) ( Supplementary Table 1 ).

To comprehend the structure of putative GmCYPs, their intron-exon organization was predicted, which revealed that GmCYPs gene structures were diverse, with the exon number varying from 1 to 9. GmCYP89A2 was predicted to be intron-less (single exon) whereas GmCYP85A1, GmCYP90A1 and GmCYP78A5 consists of 8, 7, 6 introns and 9, 8, 7 exons respectively. GmCYP93A3, GmCYP93A2 and GmCYP93A2-like genes contained 3 exons and 2 introns, while the rest of them have 2 exons with a single intron ( Supplementary Figure 2 ). This gene structure analysis helps to determine evolutionary, structural and functional characteristics of genes. Chromosomal mapping of GmCYPs illustrated that these GmCYPs were randomly placed on 13 (out of 20) soybean chromosomes. Chromosomes 8, 9 and 11 contained 2 GmCYPs, whereas the remaining 10 chromosomes contained a single gene. GmCYPs are located on chromosome arms except for GmCYP98A2, GmCYP98A2-like and GmCYP93 implying their frequent participation in recombination process ( Supplementary Figure 3 , Supplementary Table 2 ).

Conserved motif and domain analysis is a useful method for protein function prediction (Watson et al., 2005). All the listed GmCYPs were examined and analyzed with Pfam, NCBI-CDD and Inter-pro databases for the presence of functional domains; investigated the existence of the P450 domain (conserved in cytochrome family with CDD and Pfam: PF00067, Inter pro: PS00086 id). SMART and TMHMM executed the presence of trans-membrane regions in all listed GmCYPs. Additionally, all genes also harboured Cytochrome P450 Cysteine heme iron ligand signature: IPR001128, detected by ExPasy-Prosite server. In our study, MEME online tool detected the presence of 10 possible motifs ( Supplementary Figure 4 ). Motif 1, 2, 3 and 4 appeared in all candidate GmCYPs representing conserved motifs of the P450 superfamily: heme-binding region (FXXGXRXCXG), PXRX motif, I-helix oxygen-binding domain (AGxDT) and K-helix region (EXXR). Exploration of listed GmCYPs for conserved motifs was executed by MSA and identified these four conserved motifs ( Supplementary Figure 5 ).

Cis-acting elements present in the promoter region may be crucial in controlling the transcription of genes in response to environmental stressors and plant growth. A comprehensive analysis was executed to decipher cis-element analysis existing in the promoter region of 16 GmCYPs. Several cis-acting elements, including those responsive to stress, hormones and TF binding, were identified in the promoter regions of GmCYPs, where each GmCYPs harbours more than one cis-acting element ( Figure 2 ). Analysis revealed that these 16 GmCYPs harbour 88 light-responsive elements, which mainly include Sp1, sbp-CMA1c, TCT-motif, TATA-box, GATA-motif, I-box, Box 4, Gap-box, MRE, G-box, AE-box, chs-CMA2a, chs-CMA1a, GATA-motif, Sp1, GT1-motif, ATC-motif, AE-box, ATCT-motif, 3-AF1 binding site, ABRE3a, CACGTC, TCTTAC, TTACTTAA, GATAGGA motif. However, 46 phytohormone-responsive cis-elements were predicted; 16 methyl jasmonate hormone-responsive elements (TGACG-motif and CGTCA-motif), 14 ABA-responsive elements (ABRE), 10 Gibberellin hormone-responsive elements (TATC-box, CCTTTTG, P-box, GARE-motif), 3 Auxin-responsive elements (AuxRR-core and TGA-element), and 3 Salicylic acid-responsive elements (TCA-element). The presence of 3 low temperature (LTR), 22 Anaerobic conditions- responsive (ARE), 5 defense- responsive (TC-rich repeats), 5 drought-responsive elements and 1 flavonoid biosynthesis regulating (MBSI) cis-acting elements were detected. Additionally, binding sites of TFs on the GmCYPs promoter region were uncovered using PlantTFDB tool. Here, 14 types of transcription families were predicted in 11 GmCYPs promoter regions ( Supplementary Figure 6 ) while GmCYP90A1, GmCYP89A2, GmCYP78A6-like and Gm78A4-like were deprived of any TF binding sites. AP2 family were predicted in GmCYP76C1-like, GmCYP78A5-like and GmCYP98A2-like; ERF was found on GmCYP76C1-like and GmCYP82A2; Bzip was restricted to GmCYP76C1-like; NAC in GmCYP93A3, GmCYP76C1-like and Dof in GmCYP76C1-like, GmCYP93A3, GmCYP78A5-like and GmCYP85A1.

miRNA in eukaryotes is a non-coding, short RNA comprising 21-24 nucleotides that bind to and cleaves mRNA for post-transcriptional regulation (Jones-Rhoades et al., 2006). Using G. max miRNA library, miRNA target analysis for GmCYPs was performed using the psRNATarget tool; it revealed 16 GmCYPs as putative targets of 94 miRNAs when E= 5.0 ( Supplementary Table 3 ). GmCYP78A5 was targeted by the largest numbers of miRNA, 23 followed by CYP90A1 (11 miRNA), CYP93A3 (9 miRNA). In this analysis, we can conclude that miRNAs possessed diversity in regulating GmCYPs, where miRNAs inhibit GmCYPs mainly through cleavage and some via translation.

To gain insight into regulatory functions of genes controlled by a specific plant-organ, RNA-seq expression data (RKPM Values) of listed GmCYPs were procured from soybean expression explorer for investigating GmCYPs expression into six different plant-organ i.e. seed, young leaf, nodule, root, flower, pod shell, seed ( Figure 3 ). The results exhibited that all the 15 GmCYPs were expressed at least in one plant-organ, and GmCYP98A2, GmCYP98A2-like, GmCYP78A3, GmCYP78A6-like showed high expression in all six plant-organs indicating their predicted roles in soybean growth, development and countering various environmental stresses. Furthermore, some GmCYPs showed organ-specific expression patterns, such as CYP78A5 and CYP93A3, which exhibited augmented expression within the root, implying their role in the development of roots, interaction with plant microbes and stress tolerance. Similarly, CYP71AU50 and CYP85A1 were highly expressed in flowers, while CYP89A2 showed substantially greater expression in pod shells, implying their functions in reproduction/floral-related stimuli. Additionally, CYP71D10 and CYP76C1-like were found to be explicitly expressed in seeds, indicating their role in seed development and maturation. On the other hand, CYP78A4-like was expressed only in two plant-organs, seed and root, with a higher level of expression in the seed. Finally, CYP71A1 showed minimal expression in the flower and root, suggesting its possible involvement in other plant organs or stages of soybean development. However, it should be noted that expression data of GmCYP78A5 was unavailable, which limits our understanding of its organ-specific expression pattern.

Protein-protein interaction (PPI) analysis of each GmCYPs with other G. max proteins was predicted using the STRING database. PPI analysis is a comprehensive method of knowing predicted protein function. Out of 16 GmCYPs, only 8 GmCYPs have shown interactions with other G. max proteins involved in plant stress response ( Figure 4 ). GmCYP98A2 exhibited interaction with GLYMA07G02460.2 (Shikimate o-hydroxycinnamoyltransferase isoform x3) GLYMA08G23560.3 (Shikimate o-hydroxycinnamoyltransferase) and GLYMA03G27740.1 (5-O-(4-coumaroyl)-D-quinate 3’-monooxygenase); involved in phenylpropanoid pathway. GmCYP82A2 showed its association with CYP74A1 (Allene oxide synthase), which initiates the biosynthesis of jasmonic acid from hydroperoxides of free fatty acids derived from lipoxygenase, representing the first step in the process (Farmer and Goossens, 2019). GLYMA15G15830.1 (probable histone-arginine Methyltransferase 1), a positive regulator of oxidative stress tolerance, was predicted to interact with GmCYP82A2. STRING predicted GmCYP78A3 interaction with GLYMA07G02380.1 (AP2-like ethylene-responsive transcription factor), GLYMA07G20070.1 and GLYMA20G00930.1 (Checkpoint serine/threonine-protein kinase). GmCYP71AU50 exhibited interaction with GLYMA08G13230.1 (Pathogen-inducible salicylic acid glucosyltransferase), GLYMA16G02410.1 (CYP b5-like steroid binding domain-containing protein), and GPPS (Geranylgeranyl diphosphate synthase, type ii) (Barja and Rodriguez-Conception, 2021). GGPP is responsible for the biosynthesis of chlorophylls, carotenoids, or plastoquinones, which is vital for photosynthesis.

GmCYP90A1 exhibited interaction with brassinoide synthesizing enzymes; GLYMA11G07780.1 and GLYMA02G05780.2 (3-epi-6-deoxocathasterone 23-monooxygenase), GLYMA07G17410.1 (Steroid 5-alpha-reductase det2), GLYMA06G24540.1 (member of CYP450 family). It also interacts with GLYMA07G04470.1 (Typhasterol/6-deoxotyphasterol 2alpha-hydroxylase; belongs to the cytochrome P450 family) responsible for tyrosine and secondary metabolite synthesis. CYP78A5 interact with GLYMA13G07800.1 Delta(7)-sterol-C5(6)-desaturase isoform X1 and Shikimate o-hydroxycinnamoyl transferase involved in the biosynthesis of sitosterol and campesterol (Hoffmann et al., 2004). CYP78A5-like interacts with GLYMA11G37040.2 (Alpha-ketoglutarate-dependent dioxygenase alkb), GLYMA18G00970.1 (Alpha-ketoglutarate-dependant dioxygenase alkb, U) engaged in various reactions in plant metabolism such as flavonoid biosynthesis and ethylene biosynthesis (Zhang et al., 2004; Cheng et al., 2014). GLYMA16G08920.3 (heavy metal-associated isoprenylated plant protein 4) is involved in heavy metal detoxification and homeostasis mechanism. CYP85A1 interact with steroid 5-alpha-reductase (GLYMA07G17410.1 and GLYMA11G01210.1) involved in steroid biosynthesis and brassinolide biosynthesis (GLYMA02G05780.2 and GLYMA07G04470.1) and GLYMA11G26500.1 (long-chain fatty acid omega-monooxygenase), which is involved in the biosynthesis of complex polymer, sporopollenin.

Gene expression analysis of 16 GmCYPs was executed by qPCR using gene-specific primers, the results indicated that 7 of GmCYPs were upregulated while 4 of the GmCYPs were downregulated on S. litura infestation. GmCYP78A5 exhibited a maximum change in expression while Gm78A5-like, Gm98A2-like, GmCYP78A3, GmCYP76C1-like, GmCYP89A2 and GmCYP71D10 showed upregulation whereas GmCYP78A4-like, GmCYP71A1, GmCYP93A3 and GmCYP90A1 showed downregulation during infestation. However, 9 out of 16 GmCYPs were upregulated on mechanical wounding, and 3 GmCYPs were downregulated. GmCYP78A5, GmCYP89A2 and GmCYP78A5-like showed high gene expression on mechanical wounding ( Figure 5 ).

Gene expression analysis of 16 GmCYPs was performed upon the foliar spray of ethylene, salicylic acid and methyl-jasmonate to acquire insights into the role of signalling compounds in regulating the expression of GmCYPs in plant defense. On ethylene treatment, 11 GmCYPs (GmCYP82A2, GmCYP98A2-like, GmCYP78A5, GmCYP78A3, GmCYP93A3, GmCYP76C1-like, GmCYP90A1, GmCYP78A6-like, GmCYP85A1, GmCYP89A2, and GmCYP71D10) whereas 5 GmCYPs (GmCYP98A2, GmCYP78A4-like, GmCYP71A1, GmCYP78A5-like and GmCYP71AU50) were upregulated and downregulated respectively. Similarly, methyl-jasmonate application induced expression of 10 GmCYPs (GmCYP78A5-like, GmCYP98A2-like, GmCYP78A5, GmCYP78A3, GmCYP78A5-like, GmCYP71AU50, GmCYP93A3, GmCYP76C1-like, GmCYP89A2, and GmCYP71D10). However, on SA application, 4 GmCYPs (GmCYP71D10, GmCYP89A2, GmCYP78A6-like and GmCYP98A2) were induced, while 6 GmCYPs (GmCYP78A4-like, GmCYP98A2-like, GmCYP78A5-like, GmCYP76C1-like, GmCYP90A1 and GmCYP85A1) were repressed, designating the involvement of SA directly or indirectly in gene regulation. ( Figure 6 ).

Methyl-jasmonate foliar application upregulated four of the S. litura -induced GmCYPs, illustrating jasmonates’ role in regulating the GmCYPs expression upon S. litura -infestation. The presence of methyl-jasmonate sensitive cis-acting regulatory elements in the promoter region in these genes supported the notion that jasmonate is responsible for their regulation. Four GmCYPs, which had their expression elevated by ethylene treatment, were also induced by jasmonate. Jasmonate and ethylene are reported to act synergistically (Kessler and Baldwin, 2002).

Three-dimensional (3-D) structures of differentially regulated S. litura-infested 16 GmCYPs were modelled using i-TASSER, and structures were corroborated by ensuring that >90% of the protein region fell within the permitted area of the Ramachandran plot. The template, coverage (>85%), Z-score, and estimated Tm; all these parameters are listed in Supplementary Table 4 . All 16 GmCYPs chosen for analysis revealed the most typical arrangement of α-helix and β-sheets. The evolutionary conserved heme-binding region, PXRX motif, I-helix oxygen-binding domain and K-helix motifs of CYP450 proteins were highlighted using distinct colours ( Figure 7 ).

To study the enzymatic mechanism of GmCYPs, substrates of GmCYPs were retrieved from the literature survey and KEGG database. The interaction of 10 GmCYP with their specific substrate was examined through enzyme-substrate docking. The amino acid residues and polar bonds formed by substrate with GmCYPs catalyst were highlighted with magenta and yellow colors, respectively ( Figure 8 ). The global energy score of GmCYPs and substrate interaction vary from -20 to -50, indicating their strong interaction ( Supplementary Table 5 ). Earlier studies have reported that heme binding with CYP450s is essential for their enzymatic activities (Li et al., 2008). Therefore, molecular docking of a selected GmCYPs was performed with a heme ligand, and the interaction results were compared with AtCYP94C1 ( Figure 9 ). The global energy score for GmCYP98A2 like protein and heme, was -3.392, while for AtCYP94C1 and heme was -2.341.