The study involved retrieval of STEAP2-related protein sequences from the protein data bank from a protein BLAST search. The selected sequences were used for multiple sequence alignment (MSA) and phylogenetics to determine the most related species to the human STEAP2 protein. For molecular docking, a set of FDA-approved drugs were retrieved from the drug bank and used as ligands for protein–ligand docking. The protein (STEAP2) 3-D structure was first modeled and evaluated prior to the docking exercise. The protein–ligand interactions were characterized using ligplots (Figure 1).

To determine the expression of STEAP2 protein in different tissues, an expression analysis was conducted using the Human Protein Atlas (https://www.proteinatlas.org/ENSG00000157214-STEAP2) (Uhlén et al., 2015). The search term “human STEAP2” was used to find the expression of the protein in different tissues within the body.

The amino acid sequence of the query protein H. sapiens STEAP2 (Uniprot accession number: Q8NFT2) was obtained in a FASTA format from UniProt (Bateman et al., 2017). This sequence was used in the protein BLAST search for homologous sequences in the PDB.

The MSA was performed in JalView using 30 selected FASTA files from a non-redundant protein BLAST (BLASTp) search of the reference proteins in the PDB. In the BLASTp search, the algorithm parameter of maximum target sequences to return was adjusted to 5,000 and the rest of the parameters were left as default. A total of 30 sequences of the 5,000 returned were selected for MSA based on the percentage similarity in bins of 0–40, 41–60, 61–80, and 81%–100%. A number of sequences were selected from each range (Table 1).

From the table, the 0–40 and 41–60 similarity ranges took an extra hit to ensure that more distant hits from the query sequence were considered for confidence purposes.

To generate an evolution tree/dendrogram, a new alignment of the amino acid sequences from FASTA files of STEAP proteins of 30 species from the BLASTp search and subsequent 15 (only STEAP2 containing) was created via the align tab in MEGAX (Kumar et al., 2018). The sequences were copied, pasted, and alignment was performed on all using Muscle with all alignment parameters left as default. The output alignment file was used as an input for phylogenetic analysis. The maximum likelihood tree was used. In the analysis preferences, the bootstrap method was used and the number of replications were adjusted to 1,000. The (Lo and Gascuel)+discrete gamma distribution of rates across sites (G)+evolutionarily invariable sites (I) model was used and the rest of the parameter left as default. A dendrogram showing a phylogenetic relationship of STEAP2 was generated.

With the best model selected from 6 STEAP2 homology models after evaluation, the homology model by SWISS-MODEL was used for the flexible ligand, rigid receptor protein–ligand docking using a group of approved drug compounds downloaded from the drug bank. These compounds were downloaded as a SANCDB file of 2466.

The ligand and receptor molecules were prepared using prepare_ligand4.py and prepare_receptor4.py, respectively, by command line on Linux terminal. The two python scripts add polar hydrogens (Q) and define atom types (T) in the ligand and receptor files (in a PDB format), respectively, to generate their PDBQT files. The search box in the configuration file was set to dimensions x, y, and z (105 by 105 by 105 Å, with a center at: x = 122.224 Å, y = 118.600 Å, and z = 102.658 Å). These were determined using MGL tool graphical user interface with cross-sectional measurements along each axis. The median point of the line between the two furthest points of the molecule on that axis guided toward the center point to center the search box so as to cover the entire protein in all 3 dimensions.

The results of molecular docking were characterized by 3D proximity to a central residue termed as “Euclidian distance”— which is the length of a line between two points in geometrical space which can be 2 or 3D and binding energy of the ligand to the receptor molecule. All results were screened and a binding energy threshold of less than or equal to −7 and Euclidian distance of 50 Å were set. The complexes that met the selection criteria were then considered for further analysis.

The entire protein was covered by the search box during molecular docking, a method also known as blind docking to mitigate bias to any region of the receptor.

To determine the nature of the bonds formed between STEAP2 and the different ligands, Ligplot which is a program used to generate 2D representations of the nature of bonds formed between protein receptor residues and ligand atoms was used for this task (Wallace et al., 1995). The program showed the hydrogen bond interactions in green dotted lines with their respective lengths in angstroms and the different interacting protein residues in an arc shape with spokes radiating toward the atoms they are in contact with. The interactions not only accounted for hydrogen bonds but all hydrophobic interactions between the two molecules (Wallace et al., 1995).

The order of binding energy was validated using PatchDock server and AutoDock4. (Huey and Morris, 2005; Schneidman-Duhovny et al., 2005). The receptor PDB file and ligand pdbqt files were uploaded on a PactchDock server. In AutoDock4, the protocol in the tutorial by Huey et al. was followed (Huey and Morris, 2005). The binding site of STEAP2 was validated using Protein Binding Site (ProBiS) and PrankWeb webservers (Konc and Janežič, 2010a; Jendele et al., 2019). Chain A of the homo-3-mer without its HETATOMs and water molecules was the input file and was uploaded in a PDB format. All other parameters were left as default and the job was submitted for analysis at both submission webpages. The druggability of the protein was assessed using CAVIAR in the Linux terminal to identify potential cavities for ligand binding based on scores and size and ligandability among other criteria (Marchand et al., 2021). The protein chain A PDB file was used as the input file.

Messenger ribonucleic acid (mRNA) expression analysis of the FANTOM5 dataset showed that STEAP2 was most abundantly expressed in prostate tissue followed by ovarian and vaginal tissue (Supplementary Figure S2). Prostate tissue had the most STEAP2 in mean protein coding transcripts per million as compared to all other listed body tissues using the HPA dataset (Supplementary Figure S3). The analysis of STEAP2 mRNA transcripts from different tissue in mean protein coding transcripts per million in comparison to 33 other listed tissues using the GTex dataset showed prostate tissue to have the highest quantity of STEAP2 (Supplementary Figure S4). The three datasets were combined by a normalized expression level and showed the same trend (Figure 2).

A dendrogram generated using Molecular Evolution Genetic Analysis (MEGA) showed that STEAP2 was conserved mostly in species under phylum Chordata (Figure 3). The values at each node indicate the percentage times of 1,000 that particular node was redrawn as a result of genetic relationship of its taxons.

The STEAP2 containing species clustered alone in one clade and also happen to be under phylum Chordata, H. sapiens falls together with its most similar species G. gorilla (Figure 3). The rest of the clades contain other STEAP proteins and one NADP from other classes such as Pisces.

A secondary phylogenetic tree of only STEAP2s further demonstrated that the mammalian species baring STEAP2 clustered together, followed by those of Reptilia, Pisces, and Vertebrata species, A. radiata with the tree having been rooted to a yeast species S. cerevisiae (Figure 4).

The MSA image depicts the conserved residues across 15 selected sequences of STEAP2 baring species and S. cerevisiae from the original 30 from the protein BLAST search. The local regions that are highly conserved are shown with respective weblogos including the Rossman fold (amino acids residues 33–38) closely associated with the NADP domain and six other trans-membrane (TM) domains (TM1 with residues 211–232, TM4 with residues 261–284, Tm5 with residues 305–326, TM6 with residues 359–379, TM12 with residues 391–412, and TM17 with residues 432–454) associated with iron binding (Figure 5). An MSA image of all the 30 sequences is shown in Supplementary Figure.S1

Template 6HCY/6HDY (human STEAP4) was the most common top hit for five of the six modeling engines with a better coverage than 2VNS identified using PRIMO (Table 2).

The homology model from SWISS-MODEL was a homo-3-mer. Homology models from the different engines (images generated in cartoon format using PyMol) (Figure 6).

The homology model from SWISS-MODEL (Figure 6A) was a homo-3-mer made up of 437 amino acids per chain. The template used was 6HCY and the homology model had the least RMSD of 0.086 when superimposed with the 6HCY template (Table 3). Additionally, on average it had the best score by the used evaluation tools. For example, the Modeller normalized DOPE score of−5.95, the Ramachandran plot having 94.4 of its residues lying within the favored region a−3.06, and−6.33 Z-score from QMEAN and ProSA, respectively. Altogether, these results show that the SWISS-MODEL homology model was most similar model to the used template (Table 3). The comparative results of the SWISS-MODEL homology model and ab initio structural folding using AlphaFold2 are shown in Supplementary Figure S5 and a respective table. The results suggest that the homology model was a better predicted structure of human STEAP2 than the structure from AlphaFold2 that was derived by ab initio protein folding, a method better suited for query sequences with no or very low sequence identity to structurally solved sequence structures.

The rest of the results are summarized in Table 3.

From Table 3, SWISS-MODEL on average had the best homology model after validation by the listed validation tools.

Figure 7 showed the ProSA evaluation of SWISS-MODEL which was the best homology model together with its Ramachandran plot on the right. It had a z-score of−6.33.

Some of the receptor residues interacting with the drug ligands from the ligplots were found to be within highly conserved motifs of the STEAP2 protein with some forming hydrogen bonds in addition to other hydrophobic interactions. The 2D interactions are shown in panel A and the respective 3D in panel B of the images. The images of the top 3 complexes, Figures 8–10, are shown while the remaining three are attached as Supplementary Figures S6–8.

In the protein–ligand complex with Drugbank 1543 (triptorelin), Ser-372 with the hydrogen bond lies in the 6th TM domain which is part of the iron-binding domain and this is shown to be highly conserved from the MSA of only STEAP2 containing species (Figure 5). It together with many other shown protein amino acids binds with triptorelin a drug approved for the treatment of advanced prostate cancer (Figure 8).

In the protein–ligand complex with Drugbank 2 (leuprolide), the interacting protein amino acids such as Ser-372 with leuprolide lie within the TM domains that are highly conserved as shown in the MSA. These are associated with the iron-binding domain and bind to the ligand leuprolide which is also approved for the treatment of advanced prostate cancer (Figure 9).

In the protein–ligand complex with Drugbank 2328 (1,2 icosapentoyl-sn-glycero-3-phosphoserine), some of the listed residues such as Thr-210 are conserved under the NADP domain, the 5th TM domain having residues Glu-390 and Lys-303. The residues together with others form bonds with ligand Drugbank 2328 which is 1,2 icosapentoyl-sn-glycero-3-phosphoserine, an approved experimental drug for Corona virus-19 disease treatment (Figure 10).

Residues from Supplementary Figures S6 that form hydrogen bonds are all conserved highly in the NADP-binding domain and the 12th TM domain which is associated with iron binding. The ligand Drugbank 2154 is etelcalcetide, a calcimimetic drug used for the treatment of hyperparathyroidism.

The mentioned residues of Supplementary Figures S7 in TM12 are associated with the iron-binding domain of STEAP2 proteins in the MSA. The ligand Drugbank138 is lymecycline, a second generation tetracycline antibacterial used for the broad-spectrum treatment of various bacterial infections.

From the images, the protein–ligand complex of ligand 1,543 had the least number of hydrogen bonds between the ligand and protein but had the most number of hydrophobic contacts within the receptor. The summation of these bonds was, therefore, responsible for the observed high-binding energy associated with the complex (Table 4). The ligand Drugbank 1423 is droxidopa, a drug used in treatment of Parkinsonism (Supplementary Figure S8).

Table 4 shows the protein–ligand complexes that met the thresholds of the binding energy.

The complex, 01_SWISS-MODEL and Drugbank 1543 had the least binding energy.

The docking validation by PatchDock produced the same order in terms of scores based on the geometric shape complementarity of solutions. The scores were ranked in the same order as the binding energy measured using AutoDock Vina. The complex of STEAP2 and Drugbank 1543 ligand had the highest score of 11,346 followed by the complex with Drugbank 2 and Drugbank 2328 with scores 11,274 and 11,224, respectively. The least ranked complex was that with Drugbank 3300 (Table5).

Table 5 shows the top 3 solution of each complex ranked in order of their scores by PatchDock. The area stands for the approximate interface area of the complex and atomic desolvation energy is measured according to Zhang et al. (2001)

The exact docking poses from AutoDock Vina could not be reproduced by PatchDock despite maintaining the same ranking order. Another docking engine, AutoDock4, could not reproduce the exact docking poses.

The results from the ProBiS analysis of STEAP2 for potential binding sites found STEAP4 as the most locally and structurally similar protein as it was also the used template for homology modeling. The most structurally conserved residues in binding site one, for small molecules, appeared red and were in the NADP-binding domain. All six identified binding ligands were NADP’s with confidences of 4.21–3.58. Binding site two, for proteins, had all amino acids unconserved as they appeared blue and confidence from 2.42 to 2.29. No binding sites were found for nucleic acids and ions (Supplementary Figure S9A).

Consistent with these findings were those from PrankWeb whose pocket 1 was HEM binding (pocket score 24.48), 2 (pocket score 10.48) was NADP binding with the most conserved residues, pockets 3 and 6 with pocket scores 6.26 and 3.09, respectively, were also within the NADP-binding domain. Pockets 4 and 5 were all within the FAD-binding domain with scores 5.86 and 3.75, respectively. Surface pockets 4 and 5 and pocket 3 were shown to be in close proximity with docked ligands, triptorelin and leuprolide (Supplementary Figure S10).

The results from CAVIAR analysis for binding sites found 7 cavities, cavity 5 with ligandability score 0.8 lays in the most conserved residues of the chain that also bind NADP, 1 and 6 lay in the HEM-binding domain that is not very conserved with ligandability scores of 0.6 and 0.8, respectively, and 2 with ligandability score 0.8 lays in the FAD-binding domain. Cavity 2/pocket 2 was close to bound ligands triptorelin and leuprolide. Cavities with ligandability score of 0.8–1 are considered easy to the ligand and of these, cavity 5 had the best score, size (109 A), and degree of conserved residues in comparison with conservation residue score results from PrankWeb. These results are detailed in a text file of cavity identification and Supplementary Figure S9.