Fungal strains

The fungal test strains were previously isolated from lignocellulose-rich environments in central Morocco, including decaying wood, solid byproducts of olive crushing, and their compost. All strains are filamentous fungi belonging to the phylum Ascomycota. A previous study (Arif et al. 2022) screened all strains for their lignocellulosic potential by conducting qualitative and quantitative assays; several other strains were further evaluated in liquid cultures containing cellulose, lignin, and olive pomace to validate their enzymatic profiles. In a subsequent work (Arif et al. 2024), one of the most promising strains from this collection was selected and cultivated under different fermentation conditions using olive stones as a local lignocellulosic substrate. Based on these results, the following nine strains with confirmed enzyme activity and the availability of lignocellulase sequences in the NCBI database were selected for in silico analysis: A. arborescens (OPP 6), P. variotii (OPP 1), T. griseum (DWB 3), P. chrysogenum (C 9), P. crustosum (DWE 4), A. flavus (C 13), F. solani (DWI 2), P. citrinum (DWB 2), and F. equiseti (STT 4).

Retrieval of lignocellulase sequences and phylogeny construction

The sequences of cellulases and ligninases of all fungal species were retrieved from the NCBI protein database [http://www.ncbi.nlm.nih.gov]. “Fungus species” AND (“cellulase” OR “laccase” OR “lignin peroxidase” OR “manganese peroxidase”) was used as the search query. The results beyond the study’s scope were manually reviewed and filtered from the protein sequences. For each fungus species, a maximum of five sequences were randomly selected from the final list, ensuring that each sequence represents a distinct class of enzymes. Sequences were aligned using Clustal (Larkin et al. 2007) and manually filtered to exclude irrelevant results. The collected protein sequences of each enzyme group were aligned using ClustalW, and a phylogenetic tree was generated using the Neighbor-Joining method. Phylogenetic analysis was then conducted with MEGA 11 software.

Physicochemical characterization of test enzymes and secondary structure prediction

The physicochemical properties of fungal cellulases and ligninases were analyzed in silico using the ProtParam tool available on the ExPASy server (Gasteiger et al. 2003, 2005). The calculated parameters included molecular weight (MW); extinction coefficient (EC), which is important for assessing the light absorption of proteins during purification processes (Gill and Hippel 1989); theoretical isoelectric point (pI); in vivo half-life (HL), which indicates the degradation time of synthesized proteins within cells; instability index (Ii) (Guruprasad et al. 1990); aliphatic index (Ai) (Ikai 1980); Grand Average of Hydropathy (GRAVY) (Kyte and Doolittle 1982); number of negative residues (–R); and number of positive residues (+R). Additionally, secondary structural elements were predicted using the SOPMA tool (Combet et al. 2000), which enabled the analysis of secondary conformation across different segments of the protein sequences.

Homology modeling of 3D structure and model validation

The SWISS-MODEL server in the automated mode was used to predict the 3D structures of fungal enzymes (Schwede et al. 2003). The quality of the predicted models was validated according to the QMEAN score (Benkert et al. 2011). This process enabled the selection of the most appropriate enzyme model based on optimal identity and minimal energy interactions. The proposed models were further validated by uploading them to PROCHECK (Laskowski et al. 1993), ERRAT (Colovos and Yeates 1993), and the ProSA servers (Wiederstein and Sippl 2007).

Molecular docking of enzyme and substrate

Molecular docking with Autodock was performed to investigate the interaction between cellobiose and four β-glucosidases. The 3D structures of β-glucosidases from A. arborescens RYO56859.1 (PDB ID: 3zyz), P. variotii GAD99276.1 (PDB ID: 4iic), F. solani XP_046126616.1 (PDB ID: 3ahy), and F. equiseti KAJ4141610.1 (PDB ID: 5nbs) were retrieved from the Protein Data Bank (PDB) and prepared following standard docking procedures. Cellobiose, obtained from PubChem (CID: 10712), was also optimized and prepared. The docking grid boxes were centered on the predicted catalytic sites of each β-glucosidase and sized to include the main residues involved in substrate recognition. Grid box centers were set as follows: 3zyz (X = 7.3511, Y = –5.2022, Z = 13.2815), 4iic (X = 40.2146, Y = –19.5251, Z = 28.1608), 3ahy (X = –7.6994, Y = –83.1419, Z = 19.6435), and 5nbs (X = 37.3259, Y = 130.1375, Z = –5.8881). To ensure the robustness of the docking results, eight independent docking runs were performed for each enzyme-substrate system using different initial random seeds. Docking solutions were evaluated based on pose convergence, binding energy, and interaction geometry. The docking protocol and comparative analysis were further validated by including β-glucosidase B (BglB) from Paenibacillus polymyxa (PDB ID: 2O9R) as a reference structure (Khairudin and Mazlan 2013; Mazlan and Khairudin 2014). Superimposition and RMSD analysis of the predicted poses were performed using Discovery Studio Visualizer to assess pose reproducibility and binding site consistency.

Lignocellulase sequences and phylogenetic analyses

After eliminating putative, unknown, and unrelated enzymes and performing rigorous manual verification, 80 sequences were retrieved for detailed analysis (Table S1): 60 cellulases, mainly endoglucanases and β-glucosidases, and 20 ligninases, primarily laccases. No lignin peroxidase or manganese peroxidase sequences were detected in any of the test strains, and no lignocellulolytic enzymes were identified for Penicillium simile (DWE 3), Aspergillus versicolor (OPO 4), Circinella muscae (C 6), Fusarium nygamai (DWI 1), Trichoderma capillare (STT 5), Aspergillus calidoustus (OPO 6), Alternaria citri (DWE 5), Doratomyces sp. (C 10), and Penicillium sizovae (STA 3). Phylogenetic analysis was conducted on a representative subset of the sequences retrieved from the NCBI database. Distinct phylogenetic trees were generated for endoglucanases, β-glucosidases, and laccases by selecting a limited number of sequences from each species that represented the main classes of enzymes being studied.

The phylogenetic analysis of endoglucanase sequences (Figure 1A) from 7 fungal species, namely T. griseum, F. equiseti, F. solani, A. arborescens, P. chrysogenum, P. variotii, and A. flavus, revealed significant evolutionary divergence among the enzymes. The generated tree showed a well-supported clade between T. griseum and F. equiseti (bootstrap value: 100%), indicating a close evolutionary relationship between their endoglucanases. Although F. solani belongs to the same genus as F. equiseti, it occupied a separate branch, highlighting functional divergence. Similarly, species such as P. chrysogenum and A. arborescens were positioned on distinct branches, reflecting significant variations in their endoglucanase sequences. P. variotii and A. flavus, which occupied the last two branches in the tree, also showed distinct evolutionary paths, indicating differences in their endoglucanase sequences compared to those of the other species.

Figure 1

Phylogenetic analysis of sequences of (A) endoglucanases, (B) β-glucosidases, and (C) laccases of lignocellulolytic fungi retrieved from the NCBI database

As illustrated in Figure 1B, the phylogenetic tree for β-glucosidase enzymes displayed the evolutionary relationships among five fungal species: A. arborescens, P. variotii, A. flavus, Fusarium sp. (representing F. equiseti and F. solani), and P. chrysogenum. The tree revealed a well-supported clade between A. arborescens and P. variotii, with a bootstrap value of 100%. This finding suggests a close evolutionary relationship between their β-glucosidases, indicating that they may share functional similarities in lignocellulose degradation. Although A. flavus branched off separately, it had strong bootstrap support (100%), thus showing its distinct evolutionary path from the other two species, despite being grouped closer to them in this tree. At the bottom, Fusarium sp. and P. chrysogenum individually formed distinct branches, both with high bootstrap values (100%). This separation reflects the evolutionary divergence of their β-glucosidase sequences from those of the other species. The clear separation among these species highlights the diversity in their enzyme sequences, which is a key aspect for understanding how different fungi might utilize β-glucosidase in varied environmental contexts and substrates.

The phylogenetic analysis of laccase enzymes (Figure 1C) primarily focused on a representative subset of sequences rather than on the complete set of retrieved sequences to avoid generating a divergent and cluttered tree. The phylogenetic tree for laccases revealed evolutionary relationships among four fungal species: A. arborescens, P. variotii, F. solani, and A. flavus. In this tree, A. arborescens and P. variotii formed a well-supported clade with strong bootstrap support (100%), indicating a close evolutionary relationship between their laccase sequences. This finding suggests that these two species may share functional similarities in their laccase enzymes, probably playing similar roles in lignin degradation. Both F. solani and A. flavus branched off separately from the Alternaria-Paecilomyces clade; however, their respective branches were still highly supported, with bootstrap values of 100%. This observation indicates clear evolutionary divergence in their laccase sequences compared to the first clade, while demonstrating high confidence in the accuracy of their positioning within the tree.

The broad divergence observed across the fungal species suggests that each strain may have evolved specialized enzymatic mechanisms for degrading lignocellulose; this ability was likely influenced by their ecological niches. This variability in enzyme evolution indicates potential functional differences in lignocellulose breakdown efficiency.

Physicochemical properties

Table S1 presents the physicochemical properties of the studied enzyme sequences, with amino acid sequence length ranging from 224 to 891 and MW ranging from 25594.92 to 96609.56 Da. Most enzyme sequences exhibited an acidic nature, except for β-glucosidase from A. arborescens RYO56859.1, which displayed a neutral property. Additionally, endoglucanase from A. flavus KAJ1709640.1 (pI 8.18) and β-glucosidase from P. citrinum XP_056501835.1 (pI 9.36) showed alkaline characteristics.

The results for the Ii indicated that the enzymes were stable in test tube experiments and had in vivo HL up to 20 h under standard conditions (Rogers et al. 1986; Tamboli et al. 2015; Buzzo et al. 2023). The latter parameter can be predicted from the Ii (Guruprasad et al. 1990). Endoglucanases from both T. griseum (AAM77714.2 and 1W2U_A sequences) and F. equiseti AAM77703.1 were very stable with an Ii value of < 20. However, the following enzyme sequences were identified as unstable (Ii ≥ 40): A0A443HK79.1 and A0A443HK11.1 from P. variotii; KZN90305.1 from P. chrysogenum; XP_056725118.1 and KAJ5394428.1 from P. crustosum; XP_056497110.1 and XP_056501835.1 from P. citrinum; and KAJ1710886.1, KAJ1705330.1, RMZ45971.1, and RAQ73286.1 from A. flavus.

The Ai of a protein is defined based on the aliphatic side chains of amino acids such as alanine, valine, isoleucine, and leucine. It is considered a positive factor that enhances the thermostability of globular proteins (Ikai, 1980). It may also be regarded as a positive factor for increasing the thermostability of proteins, particularly globular ones. As shown in Table S1, the following proteins had the highest Ai (>75%): A. arborescens β-glucosidases (XP_028505045.1, RYO56859.1, XP_028501739.1, and RYO18646.1), A. arborescens laccases (XP_028501087.1, XP_028510645.1, RYO69350.1, RYO39750.1, and RYO39394.1), P. variotii β-glucosidases (RWQ99350.1, XP_028487745.1, and GAD93537.1), P. variotii laccases (XP_028484177.1, RWQ94532.1, and GAD97601.1), P. chrysogenum endoglucanase (KZN89060.1), P. chrysogenum β-glucosidases (KZN90305.1, KZN93408.1, KZN93114.1, and KZN91501.1), P. crustosum β-glucosidases (XP_056725118.1 and KAJ5394428.1), A. flavus cellulases (RMZ37468.1, RAQ66523.1 and RAQ63436.1), A. flavus β-glucosidases (RMZ41167.1, RMZ38224.1, RMZ37309.1, RAQ79441.1, and RAQ77926.1), A. flavus laccases (KAJ1710886.1, KAJ1705330.1, RMZ45971.1, RMZ39576.1, and RAQ73286.1), F. solani β-glucosidases (KAJ4232191.1 and KAJ4211136.1), F. solani laccases (AGT80114.1, AHD26939.1, AHZ65142.1, and AHZ65140.1), P. citrinum β-glucosidases (XP_056497110.1, XP_056501835.1, XP_056501834.1, XP_056498616.1, and XP_056495445.1), and F. equiseti β-glucosidase (KAJ4133747.1).

All the studied enzymes had negative GRAVY index values, ranging from –0.619 to –0.001. This finding indicates better interactions between these enzymes and water. The results also showed that most lignocellulolytic proteins contained more negatively charged amino acids than positively charged ones.

Secondary structure prediction and homology modeling of the 3D structure

The enzyme sequences of all strains exhibited optimum secondary structures (Table S1), except for β-turn, which showed poor conformation. The following predicted secondary structures also had a low content of α-helix: A. arborescens laccases (XP_028510645.1, RYO69350.1, and RYO39394.1), P. variotii endoglucanases (XP_028485664.1, RWQ96019.1, and GAD94262.1), multicopper oxidase domain-containing proteins (XP_028484177.1 and RWQ94532.1), A. flavus endoglucanases (KAJ1715009.1 and KAJ1709640.1), A. flavus laccases (KAJ1705330.1 and RMZ45971.1), and F. solani laccases (AGT80114.1, AHD26939.1, and AHZ65141.1). Notably, for most protein sequences, the arrangement of their different parts into random coils enhances their flexibility, providing greater conformational adaptability during enzymatic turnover (Buxbaum 2007; Tamboli et al. 2015). However, these regions are more prone to mutations and represent evolutionary hotspots (Buxbaum 2007). Based on the computed parameters, only enzyme sequences with stable physicochemical patterns and robust 2D properties were selected for further in silico analysis.

The 3D modeling of the selected sequences using the SWISS-MODEL server provided critical insights. Models with a QMEAN score of less than –4.0, indicating poor quality, were excluded, while lignocellulolytic enzymes with sequence identities greater than 50% were considered (Table S2). The QMEAN score estimates the overall “degree of nativeness” of the structural features in the model. It assesses whether the model’s QMEAN score aligns with the expected value from experimental structures of similar sizes. A QMEAN score close to zero indicates good agreement between the model and experimental structures. Only one enzyme sequence from similar models was selected for validation, ultimately yielding 23 lignocellulolytic enzyme sequences. Figure 2 shows an example of a homology-modeled 3D structure obtained from the SWISS-MODEL server.

Figure 2

Example of homology modeling of the 3D structure result retrieved from the SWISS-MODEL server for β-glucosidase RMZ38224.1 from Aspergillus flavus

Model validation and substrate docking

High-quality templates were identified based on scores of residues in the most favored region (MFR) > 90%, residues in the disallowed region (DS) < 1%, and ERRAT quality factor > 90%, reflecting good structural quality (Table S3). β-Glucosidase sequences from A. arborescens RYO56859.1, P. variotii GAD99276.1, F. solani XP_046126616.1, and F. equiseti KAJ4141610.1 were therefore considered suitable for molecular docking. Figure 3 presents the corresponding validation plots obtained from the PROCHECK, ERRAT, and ProSA servers.

Figure 3

Example of graph and plots from PROCHECK, ERRAT, and ProSA servers for 4mdo.1.A and 3zyz.1.A templates

Docking of cellobiose into the four abovementioned β-glucosidases showed consistent localization of the ligand within the predicted catalytic pocket across eight independent docking runs for each enzyme. Superimposition of the top-ranked poses revealed high reproducibility, with RMSD values below 1.0 Å, supporting the robustness of the docking protocol and the stability of the predicted binding modes. The comparative docking analysis using the reference BglB from P. polymyxa showed that cellobiose binds within the canonical active site and adopts a binding orientation identical to that observed in the four fungal β-glucosidases analyzed in this study. The docking results revealed variations in interaction patterns, reflected by predicted binding energies and hydrogen bond networks (Table S4). Among the analyzed fungal enzymes, β-glucosidase from F. equiseti exhibited a favorable interaction profile with cellobiose, with a binding energy of –5.32 kcal/mol; this interaction was characterized by multiple hydrogen bonds involving Asp88, Arg94, Lys185, Gln197, Asp276, and Gln278 residues, with distances ranging from 1.95 to 2.58 Å (Figure 4D). In contrast, the β-glucosidase from P. variotii showed a weaker predicted interaction of -0.85 kcal/mol, despite forming multiple hydrogen bonds with residues such as Arg150 (2.18 Å), Phe161 (1.77 Å), and Asp633 (1.76 Å), suggesting a less optimal spatial complementarity within the binding pocket (Figure 4B). β-Glucosidases from A. arborescens and F. solani showed intermediate binding energies of -4.97 and -4.67 kcal/mol, respectively. For A. arborescens, residues Asp61 (2.08 Å), Arg67 (1.80, 2.67 Å), Arg125 (2.12 Å), Glu441 (1.69 Å), and Lys158 (2.14 Å) formed well-positioned hydrogen bonds, with additional π-σ stacking through Trp237 (3.54 Å) (Figure 4A). F. solani interacted through Glu165 (2.15 Å) and Asp227 (2.01 Å), with weaker auxiliary contacts involving Asn301 at 4.57 Å (Figure 4C). These variations suggest that substrate recognition mechanisms differ across fungal β-glucosidases, which is primarily reflected by interaction geometry and residue-level contacts.

Figure 4

Molecular docking of (A) Alternaria arborescens (RYO56859.1), (B) Paecilomyces variotii (GAD99276.1), (C) Fusarium solani (XP_046126616.1), and (D) Fusarium equiseti (KAJ4141610.1) β-glucosidases with cellobiose as the substrate

These findings reinforce the utility of molecular docking as a predictive tool to explore enzyme-substrate interactions and identify functionally important residues (Mazlan and Khairudin 2014; Dadheech et al. 2018). Although molecular docking provides valuable insights into enzyme-ligand interactions, it merely represents a simplified approximation of enzymatic function and does not comprehensively capture in vivo complexity.

In particular, docking approaches typically neglect protein dynamics, solvent effects, and post-translational modifications such as glycosylation. Glycosylation can affect protein folding pathways, conformational entropy, and overall stability and can increase active site rigidity, leading to enhanced enzyme stability and catalytic turnover. Although glycan chains on the protein surface are often exposed and may not directly participate in substrate binding, their contribution to structural rigidity and robustness can indirectly affect enzymatic performance (Ramakrishnan et al. 2023; Hao et al. 2025). Consequently, the exclusion of glycosylation in the docking models may limit the interpretation of stability and dynamics-related features, even though the scope of this work was restricted to comparative substrate recognition. Future studies should focus on investigating glycosylation effects, for instance, through molecular dynamics simulations with explicit glycan modeling. Overall, the studied fungal species exhibit promising in silico traits that may inform their prioritization for future experimental validation in achieving lignocellulose bioconversion.