Genome-Wide Mining of Selaginella moellendorffii for Hevein-like Lectins and Their Potential Molecular Mimicry with SARS-CoV-2 Spike Glycoprotein

Multidisciplinary research efforts on potential COVID-19 vaccine and therapeutic candidates have increased since the pandemic outbreak of SARS-CoV-2 in 2019. This search has become imperative due to the increasing emergences and limited widely available medicines. The presence of bioactive anti-SARS-CoV-2 molecules was examined from various plant sources. Among them is a group of proteins called lectins that can bind carbohydrate moieties. In this article, we present ten novel, chitin-specific Hevein-like lectins that were derived from Selaginella moellendorffii v1.0’s genome. The capacity of these lectin homologs to bind with the spike protein of SARS-CoV-2 was examined. Using the HDOCK server, 3D-modeled Hevein-domains were docked to the spike protein’s receptor binding domain (RBD). The Smo446851, Smo125663, and Smo99732 interacted with Asn343-located complex N-glycan and RBD residues, respectively, with binding free energies of −17.5, −13.0, and −26.5 Kcal/mol. The molecular dynamics simulation using Desmond and the normal-state analyses via torsional coordinate association for the Smo99732-RBD complex using iMODS is characterized by overall higher stability and minimum deformity than the other lectin complexes. The three lectins interacting with carbohydrates were docked against five individual mutations that frequently occur in major SARS-CoV-2 variants. These were in the spike protein’s receptor-binding motif (RBM), while Smo125663 and Smo99732 only interacted with the spike glycoprotein in a protein–protein manner. The precursors for the Hevein-like homologs underwent additional characterization, and their expressional profile in different tissues was studied. These in silico findings offered potential lectin candidates targeting key N-glycan sites crucial to the virus’s virulence and infection.


Introduction
From the beginning of the coronavirus disease 2019  caused by the pandemic outbreak of Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) in Wuhan, China, late 2019, reports on the severe gaps in vaccine and therapeutics have come in from all over the world. With about 652 million confirmed cases and more than six million documented deaths, the disease's worldwide mortality and morbidity are still rising [1]. The evolutionary dynamics of SARS-CoV-2 epidemics vary worldwide and depend on a range of variables, including host genetics, virus genetics, co-infection, environmental factors, country health infrastructures, and public attitudes and preventative strategies [2][3][4][5][6]. However, drug discovery and immunization programs at the international level continue to take precedence to control and prevent the disease. The SARS-CoV-2 is a large, positive-sense, single-stranded, RNA coronavirus. The virus produces four major structural proteins: spike (S), nucleocapsid (N), envelope (E), and membrane (M). The spike protein mediates viral infection of cells via the transmembrane zinc protein angiotensinconverting enzyme 2 (ACE2) receptor. Structurally, the S protein is a type-I membrane trimer protein constructed from a homodimer of two subunits (S1 and S2) linked by a third membrane-embedded serine 2 protease subunit [7]. The receptor-binding domain (RBD) is located at the top of the S1 subunit. The spike protein is highly glycosylated (O-and N-glycosylation); each monomer contains 22 N-glycosylation points, partially facilitating viral virulence. Furthermore, the N-glycans at the protein's surface can mediate RBD protection against neutralization antibodies [8]. The N-glycan structures are synthesized from glycan core Glc3Man9GlcNAc2 precursor, which is then processed into three types, oligomannose (2HexNAc), hybrid (3HexNAc), and complex (<3HexNAc) structures [9]. Many drugs have been tested for their ability to inhibit the spike protein-angiotensin-converting enzyme-2 (ACE-2) interactions; they exert their potency by either blocking the RBD of the S protein [10] or obstructing the ACE2 binding domain [11]. Thus, despite enormous efforts being made for suitable candidates, there is yet a significant gap remaining to be closed in this area. While many products were developed, they were limited by global availability and/or other pitfalls in design.
Lectins are carbohydrate-binding proteins with high specificity and avidity that bind carbohydrate moieties reversibly without changing their chemical structures [12]. Over the years, plant lectins have been extensively studied for their application as antimicrobial, anticancer, anti-inflammatory, antinociceptive agents, etc. [13][14][15]. Several plant lectins have been identified as SARS-CoV-2 inhibitors. They are thought to bind the N-glycan structure near the RBD, thereby preventing the S protein from binding the ACE2. These interactions may also result in conformational adjustments that facilitate antibodies' access to the target's concealed epitopes and lead to immunological neutralization [14,16,17]. Manspecific lectins belonging to the families' legume (i.e., ConA Canavalia ensiformis and LcA Lens culinaris), Jacalin (i.e., BanLec Musa acuminata), Nictaba (i.e., Nictaba Nicotiana tabacum), Ricin-B (i.e., IraB Iris hollandica), and GNA (i.e., Gastrodianin Gastrodia elata) were reported as coronaviruses' spike protein blockers that interact with the oligo-mannose structures located near the RBD [18]. Instead of the typical carbohydrate-protein interactions, lectins can also engage in protein-protein interactions with the spike protein. A chitin-specific lectin Urtica dioica Agglutinin (UDA) isolated from the rhizomes of the Urtica dioica inhibits SARS-CoV-2 infection by binding to the RBD of the spike protein in a protein-protein manner [19].
In this in silico study, and in order to discover an inhibitory protein candidate for the SARS-CoV-2 spike protein, we analyzed members of the Hevein-like lectins (chitin-specific lectins) found in the genome of the spike moss plant (Selaginella moellendorffii, Family: Selaginellaceae), for the possibility of interacting with the spike protein's N-glycans and functioning as potential SARS-CoV-2 inhibitors. To acquire a better understanding of the physical underpinnings of the complexes, the modeled proteins were docked against the RBD of the spike protein. The complexes were then further examined using molecular dynamics simulation on Desmond and the iMODS server to provide an overview of the physical bases of the complexes. Because of its bioactive metabolites (primary and secondary), which are employed in both conventional and modern medicine as antifungal, antibacterial, antiviral, anticancer, anti-inflammatory, anti-allergy, and antioxidant agents, the spike moss was selected as a prospective source for lectin [20]. A comprehensive overview of the domain architecture, gene architectures, genomic expansion, evolutionary relationship, and expressional profiles in various tissues and organs will also be provided by the study.

Analysis of Hevein-like Gene Expansion and Evolutionary Relationship
The Hevein-like putative lectin gene IDs were checked against the Pant Duplicate Gene Database PlantDGD (available at (http://pdgd.njau.edu.cn:8080/), accessed 19 December 2022) [25] to identify the type of duplication and the duplicate genes, then the synonymous substitution (Ks) was calculated using the Ka/Ks calculator found in the TBtools v1.0986853 [26]; values higher than 1 were excluded. The phylogenetic trees were constructed using aligned trimmed lectin-like domain sequences. The maximum likelihood method and the JTT matrix-based substitution model were used in MEGA X to analyze the evolutionary relationship [27]. The final bootstrap for consensus trees was inferred from 1000 replicates.

Expressional Profile of Hevein-like Genes Based on Publicly Available Resources
The dataset of expression profiling of Selaginella moellendorffii analyzed by high throughput sequencing and deposited in the NCBI Gene Expression Omnibus (GEO) (GSE123120) was downloaded [28]. The GSE123120 file containing RNA-seq expression raw read counts per gene for semplice replicates (seeds, rhizophore, leaf, and shoot) was selected for Hevein-like genes read counts. Data were normalized, and the differential expression was analyzed by the Integrated Differential Expression and Pathway analysis (available at (http://bioinformatics.sdstate.edu/idep96/) accessed 15 April 2023) (iDEP, v0.96 [29]. Each Hevein-like coding sequence (CDS) was searched for miRNA targeting sites using the Plant Small RNA Target Analysis Server (psRNATarget, v2) under default setting (available at (https://www.zhaolab.org/psRNATarget/analysis?function=1), accessed 23 January 2023) [30].

Retrieval and Pzre-Processing of SARS-CoV-2 Spike Glycoprotein and the Lectin Structures
The 3-dimensional (3D) crystal structure of the SARS-CoV-2 S glycoprotein, at a 2.80 Å resolution (total structure weight 438.26 KDa), was downloaded from the PDB databank [42] with PDB ID: 6VXX [43]. The Schrödinger suite's Protein Preparation Wizard tool of Maestro software (Schrödinger Release 2021-2: Maestro) was used to pre-process and prepare the structure of the SARS-CoV-2 spike glycoprotein (PDB: 6VXX, designated as chain B) and the 10 Hevein-like structures. In this process, bond order was assigned, hydrogens were added, missing loops and side chains were fixed, and the uncapped Ntermini and C-termini were capped with ACE (N-acetyl) and NMA (N-methyl amide) groups, respectively. The sugar cofactor, N-acetylglucosamine (NAG), was retained in the protein structure. Afterward, the protein was optimized using PROPKA at a pH of 7.4 and minimized energy using the OPLS4 force field [44].

Identifying the Binding Sites of the Spike Glycoprotein and Lectins for Macromolecular and Ligand Docking
To identify the structural details of SARS-CoV-2 S glycoprotein and the 10 Hevein-like lectins, the NCBI Conserved Domain Search was used (available at (https://www.ncbi. nlm.nih.gov/Structure/cdd/wrpsb.cgi), accessed 18 January 2023) [45].

In Silico Molecular Docking
The HDOCK server (available at (http://hdock.phys.hust.edu.cn/), accessed 15 April 2023) [46] was used for in silico macromolecular docking of the lectins to N-linked glycans of SARS-CoV-2 spike glycoprotein using a hybrid algorithm approach. The HDOCK server uses an FFT-based method to perform global docking using the input receptor and ligand structure to sample acceptable binding conformations. Furthermore, the HDOCK docking tool performs rigid body docking by considering both the protein and the ligand (protein) onto grids. The binding energy of macromolecules was used to evaluate their binding mode and rank them based on their docking energies. In the selection process, lectin-RBD complexes were selected based on two factors: the docking lowest energy and, second, the interaction with N-acetylglucosamine (NAG) that contributed to the best complex selection when using the PDBsum: Structural summaries of PDB entries webserver (available at (http://www.ebi.ac.uk/pdbsum), accessed 15 April 2023) [47]. The HawkDock server ((http://cadd.zju.edu.cn/hawkdock/), accessed 15 April 2023) [48] was used to calculate the MM/GBSA (Molecular Mechanics-Generalized Born Surface Area) free binding energy for the docked complexes of lectins with SARS-CoV-2 S glycoprotein. The MM/GBSA was calculated based on the ff02 force field, the implicit solvent model, and the GBOBC1 model (interior dielectric constant = 1). The system was minimized for 5000 steps with a cut-off distance of 12 Å for van der Waals interactions (2000 cycles of steepest descent and 3000 cycles of conjugate gradient minimizations) and expressed in kcal mol-1. The docking modes generated were analyzed and visualized using Discovery Studio 2021 [49], which were subsequently exported as an image.
Docked lectins having intriguing interactions with the spike's RBD were subjected to further molecular docking studies against Angiotensin-converting enzyme 2 (ACE2) to evaluate their interactions at the RBD-ACE2 interface. The molecular structure of ACE2 was first retrieved from the protein complex PDB ID 6LZG (resolution 2.5) [50] and preprocessed using Maestro Software's Protein Preparation Wizard tool before being subjected to docking simulations. Protein-protein interactions (PPI) between the SARS-CoV-2 spike glycoprotein (PDB ID 6VXX) and ACE2 were investigated. This was followed by additional multiple-chain docking simulations between the RBD-lectin complex and ACE2 were performed utilizing the HDOCK server [46]. The High Ambiguity Driven Protein-protein Docking (HADDOCK) 2.4 tool was then used to perform simultaneous docking of multibody complexes between spike's RBD, lectin, and ACE2 [51]. The inputs' 3D structures were uploaded to the HADDOCK 2.4 server. The protein structural calculations were guided by defining active residues based on the previous reports for 6VXX RBD [43], Smo-lectins (Smo446851, Smo125663, and Smo99732), and ACE2 [50]. The standard HADDOCK docking protocol consists of three stages: it0, it1, and itw (rigid-body docking, semiflexible refinement in torsion-angle space, and a final explicit solvent refinement, in that order). The produced models were refined using one of three methods: energy reduction (the default in HADDOCK2.4), clustering, and ranking the HADDOCK score. The BIOVIA Discovery Studio software release 2021 [49] was used to select and analyze the 3D docked models.

In Silico Mutant Spike Protein Interactions
Five mutation points (i.e., Lys417Asn, Leu452Arg, Thr478Lys, Glu484Lys, and Asn501Tyr) are reported in spike protein's ACE 2-receptor binding motif (RBM) from many SARS-CoV-2 variants and were affected in the parental structure (PDB 6VXX) separately. The UCSF-chimera (v. 1.16) [52] command line was used to create the mutations. The mutated protein structures were then subjected to optimization and energy minimization using the protein preparation wizard tool of Maestro software (Schrödinger Release 2021-2022: Maestro). Five mutant structures were obtained, and each was docked against Smo446851, Smo125663, and Smo99732. The MM/GBSA free binding energy for each complex was calculated. Different combinations of these mutations also occur within the RBM from the variants (alpha, beta, gamma, delta, and omicron), were also used to modify the spike protein parental structure, and were docked with the mentioned lectins.

Hotspot Analysis of the SARS-CoV-2 RBD-Lectin Complex
The energetically significant hotspot residues at the interface region of SARS-CoV-2 S glycoprotein and lectins were identified using the KFC Server (Knowledge-based FADE and Contacts) (available at (https://mitchell-web.ornl.gov/KFC_Server/index.php), accessed on 27 February 2023) [53].

Normal-State Analyses via Torsional Coordinate-Association
The iMODS online server (available at (http://www.imods.chaconlab.org/), accessed on 24 March 2023) [54] was used to investigate the collective flexibility/motion functions of the lectin-spike glycoprotein complex based on normal-state analyses of their respective internal dihedral angle (torsional) coordinates [55]. Within the PDB file of lectin-spike glycoprotein docked complexes, the atoms/residues were continuously indexed. The complexes were uploaded to iMODS online server, parameters were kept as default, and the collective motions of proteins were investigated using normal mode analysis.

Molecular Dynamics Simulation (MDS)
The molecular dynamics simulation (MDS) of the docked spike's RBD with potential lectins (Smo446851, Smo125663, and Smo99732) has been performed to analyze the complex stability in the physiological environment at (100 ns) using Desmond, a package of Schrödinger Release 2021-2022 [56]. The complexes were pre-processed using Protein Preparation Wizard of Maestro, which also included optimization and energy minimization of the complexes. The systems were subjected to solvation and ionization using the System Builder tool and solvated using a TIP3P (Transferable Intermolecular Interaction Potential 3 Points) water model in a 10 Å orthorhombic box. The models were made neutral by adding counter ions where needed. To mimic the physiological conditions, 0.15 M salt (NaCl) was added. The NPT ensemble (isothermal-isobaric: moles (N), pressure (P), and temperature (T) with 300 K temperature and 1 atmospheric pressure was selected for complete simulation. The models were relaxed before the simulation. Energy minimization of the system was performed using the steepest descent method, and the OPLS_2005 force field was used in the simulation [57]. The complex's stability was assessed using 1000 frames generated from 100 ns MD simulation trajectory data. The trajectories were analyzed by plotting root mean square deviation (RMSD), root mean square fluctuation (RMSF) of the α-carbons of systems during the simulations, a radius of gyration (Rg), solvent-accessible surface area (SASA), and intermolecular hydrogen bond profile of active site of complexes. One of the main characteristics of the Hevein domains is the presence of the eight highly conserved cysteine residues. At least one domain in the Hololectins contains all eight residues (Table 1, Figure 1).   Urtica dioica agglutinin (1ENM) was used as a reference for alignment and phylogenetic analysis.

Expression Profile of Hevein-like Genes in Different Organs
The transcription profile of Hevein-like members of the S. moellendorffii varies among members within and between the same tissue or organs analyzed. For instance, the expression level of Smo125663 in seeds, leaves, and shoots is~four folds higher than in rhizophore. However, Smo99416, Smo99732, and Smo139127 expression levels are generally lower in all tissues or organs compared to the rest of the other Hevein-like genes. Moreover, the Smo99416 transcript is the only one that has a targeting cleavage site for miRNA silencing, and it is targeted by a single miRNA from the family MIR1082 (smo-niR1082a, located in scaffold 19) (see Supplementary File Figure S1 for more details).

Structural Model Building of the Lectins and Their Secondary Structures
The 3D structural models of Smo446851 and Smo443112 lectin sequences were built with the restraints from both deep learning and homologous templates of known X-ray structures of the class-1 chitinase (Glyco_hydro_19) lectin; 2DKV [58], 3W3E [59], 4TX7 [60], 6LNR [61], and 1DXJ [62] using trRosetta webserver. Moreover, trRosetta proposed other templates for Smo437354 with confidence scores > 99, but rather low identity scores in relation to the templates Expansin-like proteins (PDB ID: 3D30) [63], Beta-expansin 1a (EXPB1) (PDB ID: 2HCZ_X) [64], Pollen allergen Phl p 1 (PDB ID: 1N10) [65], and cellulose binding proteins (PDB ID: 4JS7 and 4JJO) [66]. All templates were detected by running HHsearch against the PDB70 database, and confidence, coverage, identity, E-value, and Zscore scores were reported (Figure 2, also see Supplementary File Table S1). The other seven lectin model structures were constructed based on de novo folding guided by deep learning restraints using the trRosetta webserver. The TM scores for Smo446851 and Smo443112 were 0.88 and 0.86, respectively, whereas slightly lower scores ranging from 0.86 to 0.44 were predicted for the other lectins. The TM-score is a measure of confidence estimation that falls between 0 and 1, where high scores indicate a correctly predicted topology [38] (this is detailed in Supplementary File Table S2 Figures S2 and S3A,B).
The analysis of the Hevein-like homolog sequences using SOPMA demonstrated the presence of α-helices (Hh), β-turns (Tt), random coils (Cc), and extended strands (Ec) (Supplementary File Figure S4). Almost all the lectins' secondary structures exhibited a maximum frequency of random coils (Cc), except for Smo425957 and Smo403798, indicating high protein stability and flexibility [67]. Moreover, Smo425957 and Smo403798, with a higher proportion of α-helices, accounted for 54.24% and 39.52%, respectively. Alpha helices were more likely to be found in most thermophilic proteins [68].

Molecular Docking of Lectins with SARS-CoV-2 Spike Protein
Since it has been previously established that spike protein monomer has 22 N-glycosylation 2 O-linked glycan points [15,71,72] and that Hevein-like lectins are specific for Nacetylglucosamine and chitin (polymer of GlcNAc) [73], 11 points of hybrid and complex N-glycans were of interest for the docking in this study. Due to their proximity to the ACE-

Molecular Docking of Lectins with SARS-CoV-2 Spike Protein
Since it has been previously established that spike protein monomer has 22 N-glycosylation 2 O-linked glycan points [15,71,72] and that Hevein-like lectins are specific for N-acetylglucosamine and chitin (polymer of GlcNAc) [73], 11 points of hybrid and complex N-glycans were of interest for the docking in this study. Due to their proximity to the ACE-2 RBD, the residues Asn165, Asn331, and Asn343 connected to complex N-glycan were of particular interest. The ACE-2 binding pocket and its surface N-glycans were studied for protein-protein and carbohydrate-protein interactions against each lectin, respectively.
We report that only three Hevein-like lectins ( Our in silico molecular docking showed that only three lectins (Smo446851, Smo125663, and Smo99732) were able to interact with the NAG cofactors of the spike protein S1 via hydrogen bonding. In this, two hydrogen bonds were formed between Smo446851 and NAG1307 (located at Asn343), whereas a single hydrogen bond was observed between lectins Smo125663 and Smo99732, and the N-glycan NAG1321 (located at Asn165), respectively ( Table 2). These Hevein-like lectins, which interact with the spike protein in a carbohydrate-protein manner, have also several other amino acids that interact directly with the residues located in the vicinity of the RBD. Among the 16 interacting residues of Smo446851, only Gln25 of the Hevein domain formed two hydrogen bonds with the residues from the spike protein: Cys336 and Gly339. Based on the PPI generated by the PDBSum server, 12 of the 15 interacting residues from Smo125663 located in the Heveindomains (region 19-49 and 57-91) showed good interaction with the RBD residues, of which three hydrogen bonds were formed with two residues from RBD (Ser383 and Lys386). Moreover, a salt bridge interaction was also reported between Arg91 from Smo125663 and ASP389 from the RBD. A lower number of interacting residues were observed for Smo99732 (nine amino acids), and only two residues of which Phe43 and Tyr45 are in the Hevein domain. Tyr45 shared a hydrogen bond with key residue Asn165 from the RBD (a glycosylation point for complex N-glycan) (Figure 3). Furthermore, the estimated Generalized Born Model and solvent accessibility (MM/GBSA) calculations of lectins with S glycoprotein of docked complexes provided predictions for the binding energy and detailed the contributions to binding free energy per residue to assist in the analysis of binding structures in proteins by considering the electrostatic potentials (ELE), the Van der Waals potentials (VDW), the polar solvation (GB), and the nonpolar contribution to the solvation (SA) [48]. The Smo99732, which is complexed with S protein, has a lower MM/GBSA score (−26.  • Docking score and the free energy MM/GBSA are calculated in Kcal/mol. * Confidence score > 0.7, the two molecules would be very likely to bind; between 0.5 and 0.7, the two molecules would be possible to bind; and <0.5, the two molecules would be unlikely to bind.  Initially, ACE2 was docked onto the RBD of spike protein and has shown favorable interactions with the residues of RBD Tyr505, Gly502, Gly496, Asn501, Thr500, Gln498, Tyr449, Gly446, Phe456, Lys417, Gln493, Leu455, and Tyr453 with free binding energy MM/GBSA −18.57 Kcal/mol. Afterward, a multi-body docking run conducted using the HADDOCK server revealed that the lectins are embedded between the RBD and the active site N-terminal domain (NTD) of ACE2. This leads to interference with RBD-ACE2 interactions ( Figure 4A  Smo99732 (Chain C, in red), and ACE2 (Chain C, in green). The residues in the RBD (yellow s that are involved in the interaction include Lys417, Gln493, and Asn501, while the residues NTD of ACE2 include ASP30, His34, Lys353, and Tyr354.

Hotspot Analysis of RBD-Lectin Complex
Using Knowledge-based FADE and Contacts 2 (KFC-2), the hotspot residues a clusters were predicted at the interface between S protein and lectins. The values fro KFC-a model were used as their predictability outcomes that were superior to those KFC-b model because it is based on interface solvation, atomic density, and plastici tures, as indicated earlier [53]. Specifically, residues Gln25, Asp21 from Smo4 Smo125663 residues Arg91 and Gln11, and the Smo99732 Asp25 were identified as h residues and reported higher KFC-2-a scores compared to other residues (Suppleme File Figure S5). Complimentary hotspot residues Val367, Asp364, Lys386, Gln115 Thr109 were identified as hotspot residues on S glycoprotein.

Normal-State Analyses via Torsional Coordinate-Association
The torsion angle-related normal state analysis of lectins (Smo446851, Smo1 and Smo99732) docked onto the spike glycoprotein was performed using the iM server to examine their inherited stability and conformational mobility ( Figure 5). T factor is correlated to the atom's displacement around the conformational equilib Smo99732 (Chain C, in red), and ACE2 (Chain C, in green). The residues in the RBD (yellow sphere) that are involved in the interaction include Lys417, Gln493, and Asn501, while the residues in the NTD of ACE2 include ASP30, His34, Lys353, and Tyr354.

Hotspot Analysis of RBD-Lectin Complex
Using Knowledge-based FADE and Contacts 2 (KFC-2), the hotspot residues and/or clusters were predicted at the interface between S protein and lectins. The values from the KFC-a model were used as their predictability outcomes that were superior to those of the KFC-b model because it is based on interface solvation, atomic density, and plasticity features, as indicated earlier [53]. Specifically, residues Gln25, Asp21 from Smo446851, Smo125663 residues Arg91 and Gln11, and the Smo99732 Asp25 were identified as hotspot residues and reported higher KFC-2-a scores compared to other residues (Supplementary File Figure S5). Complimentary hotspot residues Val367, Asp364, Lys386, Gln115, and Thr109 were identified as hotspot residues on S glycoprotein.

Normal-State Analyses via Torsional Coordinate-Association
The torsion angle-related normal state analysis of lectins (Smo446851, Smo125663, and Smo99732) docked onto the spike glycoprotein was performed using the iMODS server to examine their inherited stability and conformational mobility ( Figure 5). The B-factor is correlated to the atom's displacement around the conformational equilibrium. The three lectin complexes, each with S-protein, have the highest flexibility and large atomic displacements around the atomic positional range (400-600) ( Figure 5A). The complex deformability index, which indicates higher individual distortions for the docked complexes, showed general steady binding characterized by minimum deformities at coordination within the range (0-1 Å), and the peaks represent the location of hinges; the Smo99732-S protein complex is more rigid than the other two complexes ( Figure 5B). Considering that the energy required to distort the complex is proportional to its eigenvalue, the lower the eigenvalue, the easier the complex is to deform [75]. The estimated eigenvalues representing the motion stiffness for each lectin-spike protein complex were 1.92 × 10 −6 , 2.12 × 10 −6 , and 1.93 × 10 −6 for Smo446851, Smo125663, and Smo99732 complexes, respectively ( Figure 5C). The eigenvalues were inversely proportional to variance, predicting the significantly higher mobility of the lectins compared to the spike protein across collective functional motion ( Figure 5D). The iMODS server provided a covariance matrix depicting uncorrelated (white), correlated (red), and anti-correlated (blue) residue pairs. The three docked complexes have high correlated residue-pair motion compared to the uncorrelated residue-pair motion ( Figure 5E). The elastic-network model describes the lectin-spike protein complex's distinct flexibility patterns. It visualizes the atom pairs connected via springs by illustrating them according to their stiffness degrees. Dark grey is usually associated with stiffer strings. The dots indicate one spring, and a grey area indicates stiffer springs ( Figure 5F). Overall, the docked complexes showed stable binding, minimum deformities, and complex rigidity.

Molecular Dynamics Simulation
The molecular dynamics simulations were performed to evaluate the stability of

Mutant Spike Protein Interaction with Lectins
The carbohydrate-interacting lectins Smo446851, Smo125663, and Smo99732 were studied for interaction with different spike protein mutants commonly occur in SARS-CoV-2 variants (alpha, beta, gamma, delta, and omicron) (Figure 7). Five mutations related to the RBM region were selected, and five mutant 6VXX structures were produced, one for each mutation. Analyzing the 15 generated docked complexes revealed that almost all the RBM regions were now accessible for the bind of Hevein-like lectins. Additionally, many forms of interactions, including the formation of hydrogen bonds, were evident. Ile468, which is one of the residues required for the binding of the ACE-2 receptor, is mainly targeted by lectins in 11 complexes related to all mutations (Supplementary File Table S4) Figure S7). The calculated MM/GBSA free binding energy for Smo446851-mutant Glu484Lys complex was higher than reported for parent RBD and other mutant complexes (−29.16 Kcal/mol). Smo125663 showed relatively poor interaction among the mutant complexes, and was compared to its interaction with parent RBD, while Smo99732 maintained the highest free binding energy score in all mutant complexes except for the Glu484Lys-RBD complex (−8.93 Kcal/mol).   (Table 3).  • Docking score and the free binding energy MM/GBSA are in Kcal/mol. * Confidence score > 0.7, the two molecules would be very likely to bind; between 0.5 and 0.7, the two molecules would be possible to bind; and <0.5, the two molecules would be unlikely to bind. Bold entries represent the hydrogen bond formed between the carbohydrate NAG and its lectin counterpart residues.

Selaginella Moellendorffii Hevein Lectins
The spike moss Selaginella moellendorffii is a member of the lycophytes, nowadays classified under one of the three surviving families (Selaginellaceae). Its genome is the smallest reported plant genome, with a size of~100 Mbp [21]. The S. moellendorffii is a genetic model being utilized as a key to understanding how vascular plants evolved and gained many beneficial traits. Evolutionary studies showed that the S. moellendorffii genome harbors gene members belonging to all documented plant lectin families except for the Agaricus bisporus agglutinin (ABA) family [76]. The functional plasticity of plant lectins allows them to participate in various endogenous biological processes as well as defense molecules against predators and pathogens, hence their wide application in medicine and pharmaceutical fields. With the emergence of high throughput genome sequencing, the availability of genomic source databases, datasets, and their translation, and tools, novel sources of bioactive peptides and proteins were easy to identify and study. Hevein-like lectins are defense molecules; their activity is mediated by the interaction with pathogenic microorganisms' surface chitin or its derivatives N-acetylglucosamine [77]. Since the start of the COVID-19 pandemic, scientists have been searching for new drugs and vaccines to combat the virus. Many plant lectins inhibit SARS-CoV-2 by targeting the spike protein either through the interaction with the complex or high-mannose N-glycans found in the vicinity of the RBD or by interacting with the RBD residues through direct protein-protein binding [17,18]. Lectins from the Hevein (chitin-binding lectins) family mined from the S. moellendorffii genome were studied for both binding modes to the RBD of the spike protein.

Hevein Lectin Interaction with SARS-CoV-2 Spike Protein's RBD
Three Hevein-like lectins (i.e., Smo35272, Smo425957, and Smo403798) bind in the RBD region of the spike protein that spans from residues 319 to 541 in a protein-protein manner. However, none of these lectins formed any interaction with the residues that constitute the ACE2 receptor binding motif (residues 438-508, Tyr505, Gly502, Gly496, Asn501, Thr500, Gln498, Tyr449, Gly446, Tyr489, Phe456, Lys417, Gln493, Leu455, and Tyr453). These lectins bound deeply within the groove located at the top of the spike protein while some residues interacted with the amino acids flanking the RBM, especially in the region 300-399. Unlike the interaction reported for the Urtica dioica Agglutinin (UDA), examined for binding RBD (6VXX), UDA bound 14 residues (i.e., Tyr505, Gly502, Gly496, Asn501, Thr500, Gln498, Tyr449, Gly446, Tyr489, Phe456, Lys417, Gln493, Leu455, and Tyr453) which are also involved in the RBD-ACE2 interaction [19]. Therefore, the binding of Three of these amino acids, Phe20, Tyr22, and Tyr29, were responsible for the binding stabilized with hydrogen bonds during most of the simulation time [78]. These residues are conserved as part of the Hevein domain active site. However, the interaction with the RBD NAG1321 and NAG1307 via hydrogen bonds involved a different set of amino acids (Smo99732-Phe43, Smo125663-Ser30, Smo446851-Ala26, and Arg29). The Asn343, along with Asn282 and Asn331, act as a shield that covers and protects the spike protein's RBD against neutralizing antibodies. They are also involved in the binding of the ACE2-RBD. Mutations or blocking of these glycosides are reported to influence viral binding and significantly reduce its infectivity [9]. Lokhande

Interaction with Mutant SARS-CoV-2 Spike Protein's RBD
Beta, gamma, delta, and omicron are classified as variants of concern (VOC) due to their increased transmissibility, disease severity, immune response impact, and drug efficacy [79]. Several mutations occur in the RBM of the spike protein's RBD region, endowing the virus with enhanced ACE2 receptor binding and/or better antibody evasion [80,81]. By neutralizing the free virions, Triticum vulgaris agglutinin (WGA) has been shown to prevent infection of SARS-CoV-2 and its variants alpha and beta [82]. In an in vitro experiment employing Vero E6 cell lines, Griffithsin (GRFT) was found to bind to the spike protein and thereby prevent SARS-CoV-2 and its variants delta and omicron from entering the cells [83]. The NTL-125, another lectin, was investigated for its high affinity binding to the RBD of several variations. Such binding occurs through the α-helical tail of the protein, which interacts with glycan moieties and increases the binding strength [84]. Considering our findings, Smo99732 appears to be a promising candidate for carbohydrate-protein interactions with wide-type spike glycoproteins, given its complex stability, low free energy compared to other S. moellendorffii Hevein-like lectins, and capacity for improved protein-protein interactions with mutant variants' RBMs.

Conclusions
This comprehensive in silico investigation provided a novel contribution among other published studies and confirmed the effectiveness of plant mannose-and GlcNAc-specific lectins against RNA viruses, including HIV, SARS-CoV, and SARS-CoV-2. The latter two are inhibited by lectin-viral spike N-linked glycoprotein interactions. In the current study, by genome-wide search, we report on ten novel chitin-specific Hevein-like lectins from Selaginella moellendorffii, three of which, Smo446851, Smo125663, and Smo99732, were able to interact favorably with the receptor binding domain (RBD) of the spike protein. The binding specificity of these lectin homologs with spike-RBD can be further investigated using in vivo systems, such as the Vero cell line. Additionally, given the current SARS-CoV-2 mutation frequency, we anticipate that S. moellendorffii lectin-based medications, specially Smo99732, might be a candidate against newly emerging variants, such as alpha, beta, gamma, and omicron that are responsible for the rapid spread of COVID-19. Lectins are easier and highly feasible candidates than molecular vaccines, albeit toxicological analysis to address potential reactions is imperative.
Supplementary Materials: The following supporting information can be downloaded at https://www. mdpi.com/article/10.3390/cimb45070372/s1. Table S1: trRosetta templates used for lectin homology modeling; Table S2: Predicted TM-score for each lectin model; Table S3: The binding site residues of the lectins; Table S4: The hydrogen bonds formed in protein-protein interactions of mutant S protein RBD and lectins; Figure S1: Expression profile of Hevein-like genes from Selaginella moellendorffii in different tissues/organs; Figure S2: The Ramachandran plots for the Hevein-like lectin model structures; Figure S3A,B: The quality and Z-score of the Hevein-like lectin structural models (five models); Figure S4: Prediction of the secondary structure of Hevein-like lectins from Selaginella moellendorffii. Hh: α-helices, Tt: β-turns, Cc: random coil, and Ec: extended strand; Figure S5: Hotspot residue scores predicted by KFC-2 for the PPI involving SARS-CoV-2 S glycoprotein; Figure  S6

Informed Consent Statement:
No specific permits were required because all repository data servers were publicly accessible. Data Availability Statement: All data generated and analyzed during this study are included in the main article, and its supplementary data are provided in the supplementary section of this manuscript. Genomic data of Selaginella moellendorffii were publicly available and freely accessible, and all related UTR links were provided within the article under relevant mentions.

Conflicts of Interest:
The authors declare no conflict of interest.