Reconstruction of Ancestral Sequences of Vg

Our team initially utilized the BLAST function in the NCBI database to identify protein sequences with high similarity to the overall sequence of Vg. We obtained 50 protein sequences from various species, all exhibiting a similarity greater than 71.3% to Vg. After removing 17 duplicate records, we aligned the remaining 33 protein sequences using the ClustalW function in MEGA X software (Figure 1 A) and constructed a corresponding phylogenetic tree using the maximum likelihood algorithm based on the Jones-Taylor-Thornton (JTT) model in MEGA X software (Figure 1 B).

We conducted an evolutionary analysis of Vg based on the phylogenetic tree and homology comparison results (Figure 1 A, B). From the perspective of natural evolution, conserved domains and mutation sites were identified (Figure 1 A). By examining different positions within the mutation intervals, we categorized the mutation sites into two types: functional domain mutation sites and non-functional domain mutation sites (Figure 1 C). Since the pharmacology of the Vg we developed is associated with the VDAC-1 receptor protein, the functional domain corresponds to the region that directly interacts with the receptor (Figure 1 D). Therefore, we focused our mutation studies solely on the functional domain mutation sites.

Figure 1. Evolution and homology analysis of Vg based on reports.(A) Utilized the align by clustalw function of megax software, where the yellow regions represent conserved domains and the white areas indicate mutation sites. (B) The data in this section is sourced from the blast in the NCBI database, employing the maximum likelihood algorithm of the Jones-Taylor-Thornton (JTT) model in megax software to construct a phylogenetic tree, which was beautified using cnsknowall (https://cnsknowall.com/). The red star marks Vg, the species in the green branches belong to Fabaceae, and those in the blue branches belong to Euphorbiaceae. No protein sequences with high overall similarity were found in animals and microorganisms. (C) The conserved domain database (CDD) from the NCBI database was used, and through CD-search, a highly reliable albumin-I superfamily functional domain of Vg was identified. According to the annotations from the conserved domain database (CDD), the albumin-I superfamily functional domain can bind to the 43 kDα receptor protein VDAC-1. (D) The image depicting the binding of Vg to the receptor VDAC-1 was generated using Autodocking Vina, showing the conformation when the free binding energy of Vg to the receptor VDAC-1 is minimized at -8.9 kcal/mol. The visualization was enhanced using pymol, where VDAC-1 is displayed in surface form, the functional domain of Vg is shown in mesh form, and the non-functional domain of Vg is represented in cartoon form. The figure clearly illustrates the direct interaction of the functional domain of Vg with the molecular binding pocket of the receptor. VDAC-1 is displayed in surface form, the functional domain of Vg is shown in mesh form, and the non-functional domain of Vg is represented in cartoon form. The figure clearly illustrates the direct interaction of the functional domain of Vg with the molecular binding pocket of the receptor.

Mutation Site Strategy of Vg Guided by Machine Learning

We plan to subsequently synthesize Vg in Pichia pastoris eukaryotic heterologous expression system, which is quite different with plant expression systems. In order to achieve successful expression of Vg in Pichia pastoris, we employed the artificial neural network (ANN) algorithm within the netphos3.1 algorithm to analyze and predict the key domain in all the selected sequences. This approach leverages extensive and systematic learning to capture the relationship between protein three-dimensional structure and function, enabling precise prediction of phosphorylation sites^[1] (Figure 2 A). By processing the results derived from the NetPhos3.1 algorithm (Figure 2 B), we identified phosphorylation sites at positions 18 and 36, which can influence the intracellular localization of the protein. Phosphorylation can alter the hydrophilicity of the protein or its interaction with the cell membrane, leading to the translocation of the protein from the cytoplasm to the cell membrane or nucleus, which is crucial for cellular function and signal transduction. We confirmed the presence of phosphorylation sites at positions 18 and 36, which can affect the protein's localization within the cell. Through changes in the protein's hydrophilicity or interactions with the cell membrane, phosphorylation can result in the movement of the protein from the cytoplasm to the cell membrane or nucleus, which is vital for cellular functions and signaling pathways^[2].

The selection of mutation sites mentioned above is based on the arrangement of amino acids, considering the primary structure level. Therefore, our team utilized the deep learning architecture and neural networks of AlphaFold3 (Figure 2 C) to predict the three-dimensional structure of Vg, employing Pymol to visualize the three-dimensional conformation of Vg. Since Vg is a polypeptide without a quaternary structure, we will select mutation sites based on the secondary and tertiary structures.

We predicted the secondary structure of Vg using the latest jpred4 algorithm from the jnet algorithm series. The results indicate that Vg has three β-sheet regions at positions 21 to 31, while other regions consist of β-turns and random coils (Figure 2 D). In the annotation of the albumin superfamily functional domain of Vg, there is a correlation between the function of this domain and the three β-sheet regions (Figure 1 C), which allows us to narrow the research area to the mutation sites from 21 to 31.

In the three-dimensional structure of Vg, the significant advantage of Vg as an oral hypoglycemic agent is primarily considered, a property determined by its globular three-dimensional structure. The key to Vg forming a stable globular structure lies in three pairs of disulfide bonds (Figure 2 E). Observations reveal that disulfide bonds exist between amino acids at positions 7 and 22, 3 and 15, as well as 20 and 32. Consequently, the sites that form disulfide bonds in the mutation sites were eliminated.

Finally, considering that the essence of biological evolution is gene mutation, we performed tblastn operations on 33 proteins with high overall similarity to Vg in the NCBI database for gene prediction. For the predicted genes, we utilized the branch-site model in the PAML software, employing likelihood ratio tests to analyze the accelerated changes in nonsynonymous substitutions in specific genes, obtaining values for the constant. The identified positive selection sites were then traced back to the proteins, serving as important references for valuable mutation sites in proteins. For the predicted genes, we utilized the branch-site model in the PAML software, employing likelihood ratio tests to analyze the accelerated changes in nonsynonymous substitutions in specific genes, obtaining values for the constant. The identified positive selection sites were then traced back to the proteins, serving as important references for valuable mutation sites in proteins (Figure 2 F)

Figure 2. Machine learning-assisted selection of Vg mutation sites.(A) Illustrates the architecture of the artificial neural network (ANN) algorithm, which simulates the cognitive and learning processes of the human brain, showcasing a large and systematic learning framework that captures the relationship between protein three-dimensional structure and function. (B) NetPhos3.1 exports predicted phosphorylation sites, with points marked as "yes" or exceeding the purple boundary indicating phosphorylation sites predicted with higher confidence by machine learning, resulting in predicted phosphorylation sites at positions 18 and 36. (C) Displays a simple architecture of AlphaFold3, divided into three modules: blue represents MSA & template, green represents PairFormer, and pink represents the diffusion model. The data transmission route is divided into three parts: yellow indicates input data, blue indicates abstract network activations, and green indicates output data. (D) The JPred4 algorithm exports predicted secondary structures of Vg, where the green regions indicate β-fold regions, with β-fold regions present at positions 21-31. The data transmission route is divided into three parts: yellow indicates input data, blue indicates abstract network activations, and green indicates output data. (E) Displays the three-dimensional structure of Vg predicted by AlphaFold3, with disulfide bonds colored red using PyMOL software. (F) Shows the branch-site algorithm of PYMAL.

Through the guidance of the aforementioned machine learning, we developed mutation strategies and plans, focusing on the evolutionary strategies of different amino acids at positions 25, 27, and 28 (Figure 3 A, B, C). We combined various mutation schemes at these sites, resulting in a total of 26 mutant products, which we named (R25A, R25A_V27F, R25A_V27F_V28A, R25A_V27F_V28I, R25A_V27L, R25A_V27L_V28A, R25A_V27L_V28I, R25A_V28A, R25A_V28I, R25G, R25G_V27F, R25G_V27F_V28A, R25G_V27F_V28I, R25G_V27L, R25G_V27L_V28A, R25G_V27L_V28I, R25G_V28A, R25G_V28I, V27F, V27F_V28A, V27F_V28I, V27L, V27L_V28A, V27L_V28I, V28A, V28I), and utilized AlphaFold3 to depict the three-dimensional structures of these various mutant products.

Figure 3. Mutations site in Vg.(A), (B), and (C) respectively display the mutation sites 2, 5, 2, 7, and 2, 8 of Vg. The sites 2, 5, 2, 7, and 2, 8 are represented in stick form, while the remaining parts are depicted in cartoon form.

Virtual Screen for Vg Binding Molecular

Our team plans to first evaluate the binding potential of molecular ligands for Vg through a machine learning-based empirical scoring mechanism, distinguishing between active and inactive compounds to eliminate some mutation schemes that become inactive due to random mechanical combinations. We utilized the Chemscore algorithm in Maestro software, which decomposes the total energy of the protein-ligand complex into multiple energy terms, removing unreasonable geometric configurations and mutual positions, ultimately obtaining the activity of the molecular ligands through a scoring function. Subsequently, our team assessed from the perspective of the stability of the protein-ligand complex, selecting the Hex scoring function from Discovery Studio software. By collecting and calculating various energy terms of the protein-ligand complex after virtual docking, we weighted and combined these terms using the scoring function. Finally, to achieve the speed of reaching the binding energy equilibrium of the protein-ligand complex, we selected the Glide module in the Glide score of Maestro software, employing the OPLS all-atom force field. We collected data derived from the three methods using the all-atom model mechanical analysis in the OPLS all-atom force field, processed the data, and plotted a three-dimensional graph using Origin (Figure 4 B).

We proceeded to use protein-ligand blind docking to simulate molecular behavior in a real biological environment. By analyzing the results of molecular docking, we identified the binding pocket of the receptor VDAC-1 (Figure 4 A), which provided an important parameter for the next step of specific docking. In the protein-ligand blind docking, we selected the Vina force field in conjunction with the AutoDock algorithm. We utilized the Vina force field to derive a mathematical expression^[3] that describes the relationship between the system's energy and its particle coordinates, calculating the total energy of the molecular system to assess the stability of different conformations. This allowed us to obtain the three-dimensional structure of the most stable conformation and the system energy, which served as an important parameter for the next step of specific docking. Finally, we employed the global optimization algorithm of AutoDock to draw conclusions. We also utilized the same MMGBSA calculation scheme along with the TIP4P water and solvent model, and the GBSA solvent model, using different molecular force fields, specifically the Amber force field and the OPLS-AA force field. The free binding energy calculation results from the two force fields were designated as Amber and Prime MMGBSA. We collected data derived from the three methods and processed it, using Origin to create a three-dimensional graph (Figure 4 C).

Figure 4. Virtual Screen for Vg binding moleculars.(A) First, the binding pocket of vdac-1 was predicted using the sitmap function in Maestro software, resulting in ten binding pockets. The binding site of vdac-1 with Vg was then determined based on the binding positions obtained from the autodocking blind docking algorithm. In the figure, vdac-1 is displayed in stick form, the binding sites are shown as white surfaces, and the binding pockets are represented as purple cubes. (B) The x-axis represents the results derived from the chemscore algorithm in Maestro software, the y-axis represents the hex scoring function from Discovery Studio software, and the z-axis represents the glide score results from Maestro software. The binding site of VDAC-1 with Vg was again determined based on the binding positions obtained from the autodocking blind docking algorithm. In the figure, vdac-1 is displayed in stick form, the binding sites are shown as white surfaces, and the binding pockets are represented as purple cubes. The red spheres indicate Vg, while the blue spheres represent the other 32 evolutionary strategies. In the data processing, since the chemscore results are within the range of (0-1), and both the hex scoring function and glide score represent binding free energy, we artificially took the negative values of all binding free energy data for easier visualization and observation of conclusions. (C) The x-axis represents the results derived from the vina force field combined with the autodocking algorithm, the y-axis represents the results from the mmgbsa calculation scheme using the opls-aa force field, and the z-axis represents the results from the mmgbsa calculation scheme using the amber force field. The red spheres indicate Vg, the orange spheres represent the top three mutation strategies in terms of binding effectiveness with vdac-1 excluding Vg, and the blue spheres represent the other 29 mutation strategies. In the data processing for figures (B) and (C), we took the negative values of all binding free energy data to maintain consistency.

After six rounds of calculations, two methods, the empirical function and mmgbsa, utilized four types of molecular force fields: the OPLS all-atom force field, the Vina force field, the Amber force field, and the OPLS-AA force field. This comprehensive evaluation demonstrated that Vg has a significant advantage in binding to the receptor VDAC-1 compared to various evolutionary strategies, allowing for the direct use of the wild-type Vg in the subsequent exploration of wet lab synthesis methods.

To correctly express Vg, a hypoglycemic polypeptide derived from soybean consisting of 37 amino acids. We chose Pichia pastoris GS115 as our chassis microorganism. Pichia pastoris is a kind of yeast that can utilize methanol as its sole carbon and energy source. It offers significant advantages over other expression systems in terms of protein processing, secretion, post-translational modification, and glycosylation^[4]. To achieve our goal, we have constructed three plasmids and then transformed them into GS115 for expression. The plasmid, parts we used are shown in parts.

pPIC9K-his-DDDDK-Vg

Our goal is to consistently produce Vg using GS115. As mentioned in the engineering page. We already constructed three plasmids in our project. For expressing Vg successfully, the vector we selected is pPIC9K, which codes a signaling peptide that makes the expressed protein secret into the culture environment. There is an AOX 1 promoter on the pPIC9K vector, encoding an alcohol oxidase that breaks down methanol as a carbon source^[5].

Figure 5. Visualization of the pPIC9K-his-DDDDK-Vg for overexpression.Recombinant plasmid profile of PIC9K-3xhis-DDDDK-Vg. The expression of P37 gene in the recombinant plasmid was strictly regulated by methanol using AOX promoter, and α factor secretion signal was an extracellular signal peptide mediating Vg fermentation and secretion into culture medium. Ek (Enterokinase) can recognize DDDDK sequence in protein efficiently and specifically, and separate P37 from other non-target proteins to obtain P37 with high purity^[6].

Since we required a large number of plasmids, we transferred the recombinant pPIC9K plasmid to DH5α to obtain a large number of plasmids. After the plasmid was extracted, the plasmid was linearized by single digestion with Sac I (Figure 6 A), and we transferred the linearized plasmid into Pichia pastoris by chemical transformation. After the dilution of the transformed bacterial solution, we cultivated the single colonies. Meanwhile, the positive clones were further selected by dual resistance to kanamycin and ampicillin. Subsequently, we first performed colony PCR to select the positive clones. We designed the PCR results for two pairs of primers (Figure 6 B), and all the selected plasmids were at the expected position, consistent with the position of the positive control.

After obtaining the validated positive clones, we induced the expression of the target protein by methanol, and used Western blotting (Figure 6 C) and Coomassie Brilliant Blue staining (Figure 6 D) to verify the successfully expressed recombinant protein in the supernatant.

Figure 6. Biosynthesis of Vg in Pichia pastoris.(A) Linearization of pPIC9K-his- DDDDK-Vg recombinant plasmid. The pPIC9K-his-DDDDK-Vg recombinant plasmid was linearized by Sac I, and it was identified that the dissected band should run slower with the intact plasmid running ahead, consistent with the theory. Lane Sac I was the experimental group of linearized pPIC9K-3xhis-DDDDK-Vg recombinant plasmid, and Vector was the group of undigested plasmid. (B) Screening positive clones stably expressing Vg. The extracted Pichia pastoris genomic DNA was analyzed with vector universal upstream and downstream primers 5, AOX-F, and 3, AOX-R about 670 bp in size (up panel) , and P37-F and 3, AOX-R, upstream and downstream of the target fragment, about 300 bp in size down panel . Except sample 10, the other samples are in accordance with expectation. 1-20 was the recombinant GS115 strain, G was the negative control of Pichia pastoris strain without transforming P37 plasmid. P was the positive control using pPIC9K-his-DDDDK-Vg as the template. (C) The expression of P37 was detected by Coomassie Brilliant Blue staining. The difference between GS115 and the culture supernatants of strains 2,3,5,6,7 was observed, and the size was inconsistent with the theory(16.9 kDα). G was negative control group of Pichia pastoris strain without transforming P37 plasmid.

After the recombinant protein was successfully obtained, we performed nickel column affinity chromatography on the supernatant of culture medium to explore whether His-tag could bind to nickel column successfully. It was found that the recombinant protein could not be attached to the nickel column. We suspect that the signal peptide at the front of the nickel column may mask the binding of His-tag to the nickel column (Figure 7).

Figure 7. Vg enrichment and purification with Ni-NTA.Proteins from culture supernatant (S), flow through (FT), wash (W) obtained by Ni-NTA affinity chromatography were separated by SDS–PAGE. M protein markers (kDα). 10 nm, 10 mm, 50 mm, 200 mm, 500 mM means that we used the corresponding concentration of imidazole to wash the column, and then washed the protein with low concentration [(10 mM) of imidazole first, repeated three times, and then eluted the target protein with high concentration gradient (50 mM, 100mM, 200mM of imidazole)].

pPIC9K-3Vg-His

To prevent His-tag from being masked, we moved it to the C-terminal and ran the experiment again. By reviewing the literature, we learned that we could increase the expression level of our Vg by increasing the copy number of the target gene, and therefore, we changed the original single-copy Vg to three copies^[7]. We also learned that, in addition to DDDDK being suitable for protein purification, Asp-Pro is also a good choice and is more economical and efficient to use^[8,9]. So our plasmid design turns out to be like this (Figure 8).

Figure 8. Visualization of the pPIC9K-3Vg-His for overexpression.We were able to send the newly designed plasmid to the company for synthesis, and this is our plasmid sequence (Figure 9 A). Then we carried out Sac I digestion (Figure 9 B). The digested plasmid was transferred into Pichia pastoris by chemical transformation, and the obtained recombinant bacterial solution was diluted and coated. The positive clones of GS115 strain 1,2,4,7,8,9 were selected and identified by PCR (Figure 9 C).

Figure 9. (A) the peak diagram of the recombinant plasmid pPIC9K-3Vg-His, which is consistent with the theoretical sequence as shown in the peak diagram. (B) linearization of pPIC9K-3Vg-His recombinant plasmid. The recombinant plasmid pPIC9K-3Vg-His was linearized at Sac l site. It was identified that the cut band should run slowly and the intact plasmid ran ahead, which was consistent with the theory. Lane Sac l was the experimental group of linearized pPIC9K-3Vg-His recombinant plasmid, and Vector was the control group without restriction enzyme digestion. (C) Screening positive clones stably expressing Vg. The extracted Pichia pastoris genomic DNA was used to judge whether Pichia pastoris expressed P37 with primers upstream and downstream of the target fragment P37-F and P37-R, and the size was about 1,900 bp (Figure 5, C, Left). Using vector universal upstream primer 5, AOX-F and target gene downstream primer P37-R, one primer carrier, one primer on the target gene, to determine whether the target gene is transferred to Pichia pastoris, about 570 bp in size (Figure 5, C, Right). 1-9 is the recombinant GS115 sample expressed by Pichia pastoris strain. G is the negative control of Pichia pastoris strain without transforming P37 plasmid.

After obtaining the validated positive clones, we induced the expression of the target protein by methanol, and used Western blotting (Figure 10 A) and Coomassie Brilliant Blue staining (Figure 10 B) to detect expressed recombinant protein in the supernatant. Unfortunately, we still have not successfully detected the target band, considering that our recombinant protein may be too small molecular weight, not easy to express or detect.

Figure 10. (A) the expression of Vg was detected by Coomassie Brilliant Blue staining. There was no significant difference between GS115 and the culture supernatant of experimental group. (B) Vg expression was detected by Western blotting. The results of SDS-PAGE electrophoresis and His antibody incubation showed that the recombinant plasmid was not expressed or detected in Pichia pastoris. His was negative control group, 1-6 was recombinant Vg sample expressed by Pichia pastoris strain, G was negative control group of Pichia pastoris strain without transforming Vg plasmid.

pPIC9K-NK-3Vg-His

We found nattokinase (NK), an alkaline serine endopeptidase secreted by Bacillus subtilis (natto), has the advantages of long half-life, high specificity, low side effects and direct oral administration with molecular weight of 27.7 kd^[10] . We added it to the original plasmid, which can compensate for the low molecular weight defect of our target protein. At the same time, in the process of late purification, we can collect Vg and NK, respectively, as our expression products (Figure 11).

Figure 11. Visualization of the pPIC9K-NK-3Vg-His for overexpression.

Then we carried out Sac I digestion (Figure 12 A). The digested plasmid was transferred into Pichia pastoris by chemical transformation, and the obtained recombinant GS115 solution was diluted and coated. The positive clones of GS115 strain 1,3,8 were selected and identified by PCR (Figure 12 B).

Figure 12. (A) pPIC9K-3Vg-His recombinant plasmid linearization. The recombinant plasmid pPIC9K-3Vg-His was linearized at Sac I site. It was identified that the cut band should run slowly, and the intact plasmid ran ahead, which was consistent with the theory. Lane Sac I was the experimental group of linearized pPIC9K-NK-3Vg-His recombinant plasmid, and Vector was the control group without restriction enzyme digestion. (C) Screening positive clones stably expressing P37. The genomic DNA of Pichia pastoris was extracted to examine the signature sequence with universal upstream primers (AOX-F and AOX-R), and the size of product was about 1700 bp (b, upper panel). The PCR product amplificated by primers Vg-F and AOX-R, was about 1300 bp in size (b, lower panel), which is indistinguishable from the expected sample 1,3,8. 1-20 was the recombinant P37 expressed by Pichia pastoris.

After obtaining the validated positive clones, we induced the expression of the target protein by methanol, and used Western blotting (Figure 13 A, B) to verify the successfully expressed recombinant protein in the supernatant. Finally, we successfully constructed the GS115 expression system of Vg. The effect of P37 concentration is obvious, so we can extract P37 from large-scale induction culture.

Figure 13. (A) Detection of P37 expression by Western blotting. SDS-PAGE electrophoresis and His-antibody incubation showed that there was a weak signal band between GS115 and the supernatants of No. 1,3,4,5,7 in 50 kDα. 1-6 was recombinant P37 expressed by Pichia pastoris strain. G was negative control group of Pichia pastoris strain without transforming P37 plasmid. (B) P37 expression was detected by Western blotting after concentration.

Conclusion

Our team has noted that the mutation at positions 28 to 'A' holds significant potential. In the future, we can further investigate these positions using synthetic biology methods to explore more theoretical mechanisms of molecular binding. In selecting mutation strategies, our team utilized the local alignment algorithm of the NCBI BLAST function. Subsequently, we can learn from the profile-HMM algorithm mentioned by the 2023 NUU-China team, which will enable us to identify more proteins that have different structures but similar functions. Through comprehensive analysis, we aim to uncover the essential mechanisms of the docking between Vg and VDAC-1.

The aim of our project is to construct a microbial expression system that can efficiently express peptide Vg. Our major achievements and findings include the successful construction of a Vg-expressing GS115 strain, which contributes to the industrialization of Vg. Our achievements can enable the world's diabetes patients to enjoy the Vg-related treatment, and improve their daily lives.