In Silico Prioritization of Plant Growth Regulators as Candidate Modulators of Microalgal Lipid Biosynthesis for Biofuel Production

1 Introduction

Microalgae, comprising prokaryotic cyanobacteria and eukaryotic photosynthetic protists, are the most abundant primary producers and play a fundamental role in aquatic food webs [1]. They contribute to nearly half of global oxygen production and primary productivity. Owing to their robust photosynthetic ability, rapid growth, and efficient nutrient conversion, microalgae are increasingly valued across diverse industries [2]. The increasing demand for sustainable energy has intensified interest in alternative biofuels, with microalgae recognized as promising candidates [3]. These organisms can convert atmospheric carbon dioxide into lipid-rich biomass suitable for biodiesel and other biofuels. Biodiesel production relies on polyunsaturated fatty acids (PUFAs) and triacylglycerols (TAGs) accumulated in microalgal cells, yet maximizing lipid yields remains technically challenging [4]. Strategies such as genetic engineering, culture optimization, and nutrient stress induction have been explored [5], but often face limitations including higher costs, impaired growth, or unstable yields [6].

Phytohormones (PH) and different plant growth regulators (PGRs), such as auxins, cytokinins, gibberellins, abscisic acid, ethylene, salicylic acid, polyamines, betaines, jasmonates, brassinosteroids, and strigolactones, have been extensively studied in higher plants for their roles in enhancing plant yields and stress responses [7]. Recent reports suggest that exogenous application of these compounds can also enhance microalgal biomass accumulation and influence metabolic pathways, promoting metabolite biosynthesis while maintaining cell viability [8]. Such treatments may affect enzymatic activities in lipid metabolism, increase photosynthetic efficiency, or trigger adaptive stress responses that favor lipid storage [2]. Physiological studies have documented increased lipid productivity in microalgae treated with, auxins (IAA, IBA), cytokinins (K, cZ, 6-BAP), gibberellins (GA3), strigolactones, abscisic acid, methyl jasmonate, brassinolide (BL), epi-BL, and salicylic acid [9].

Despite these observations, many PGRs remain underexplored, with more than 70 brassinosteroids [10] and over 250 gibberellins described. Auxins include several indole derivatives, and cytokinins encompass multiple adenine analogues. A cost-effective and rapid screening strategy is therefore needed to prioritize compounds for in vivo assays [11]. Computational modeling offers a complementary approach, enabling the prediction of candidate binding interactions and electronic properties without the resource demands of experimental trials. Such insights may guide the rational selection of PGRs for further study in lipid biosynthesis and accumulation.

This study represents an in silico exploration of sixty-five PGRs as candidate binders of microalgal lipid biosynthesis enzymes. Five representative proteins were selected–FabD, KASII, FabG, FATA, and GPAT–covering initiation, elongation, reduction, termination, and assembly steps of fatty acid and triacylglycerol pathways. An integrative computational workflow was applied, combining virtual screening, molecular docking, molecular dynamics (MD) simulations, and density functional theory (DFT) calculations. The pipeline is designed to generate hypotheses by identifying predicted binding modes, candidate interactions at catalytically relevant regions, and supplementary electronic descriptors. These computational findings are intended to prioritize ligands for subsequent enzymatic and cellular validation, with the broader aim of informing strategies to enhance microalgal lipid productivity for biofuel applications.

2 Methods

2.1 Preparation of Ligands

Ligand structures were retrieved from the PubChem database [12]. A total of sixty-five plant growth regulators representing eleven chemical classes (Table 1) were selected to ensure structural diversity. Each compound was initially pre-optimized using the MMFF94 force field. Conformational searches were performed, and the lowest-energy conformer was retained. Final geometry optimization was performed at the same level of theory as implemented in AVOGADRO software [13]. The optimized structures were subsequently used as input ligands for molecular docking simulations against the selected microalgal enzymes.

2.2 Preparation of Target Enzymes

Five enzymes were selected from the Protein Data Bank (PDB) to represent key steps in fatty acid synthesis (FAS) and glycerolipid assembly pathways. The crystal structures retrieved from the RCSB database included: malonyl-CoA:ACP transacylase (FabD; PDB: 4RR5) [14], β-ketoacyl-ACP synthase II (KASII; PDB: 1E5M) [15], 3-oxoacyl-ACP reductase (FabG; PDB: 4DMM) [16], acyl-ACP thioesterase (FATA; PDB: 9MQF) [17], and glycerol-3-phosphate acyltransferase (GPAT; PDB: 8IA1) [18]. These enzymes cover initiation, elongation, reduction, termination, and assembly steps, providing a comprehensive view of the plastidial Type II FAS system found in green algae (Chlamydomonas reinhardtii, Myrmecia incisa) and cyanobacteria (Synechocystis, Synechococcus). Structural parameters and biological origins are summarized in Table S1 (Supplementary Data). Protein structures were prepared using BIOVIA Discovery Studio Visualizer 2024 [19] by removing heteroatoms, co-crystallized ligands, and solvents to optimize docking conditions. Hydrogen atoms were added, and the charges were adjusted. Final corrections were performed using Swiss-PdbViewer [20].

2.3 Virtual Screening and Molecular Docking

Virtual screening based on molecular docking was performed to prioritize candidate ligands among the sixty-five plant growth regulators. Docking was carried out against the active site of the five target proteins involved in lipid metabolism using AutoDock Vina, as implemented in the PyRx software platform [21]. The search center (SC) coordinates and box dimensions (BD) applied in docking are provided in Table S2 (Supplementary Data). To validate the docking protocol and establish benchmarks for binding affinities, reference ligands were included alongside the screened compounds. These reference molecules comprised natural substrates (malonyl-CoA for FabD, NADPH for FabG, and sn-glycerol-3-phosphate for GPAT), substrate analogues (PN7 hexanoyl-mimetic for KASII and oleic acid for FATA), and known inhibitors (corytuberine for FabD, platensimycin for KASII, luteolin for FabG, WP2 spirolactam for FATA, and FSG67 for GPAT). All reference ligands were retrieved from their corresponding PDB entries or PubChem and prepared using the same protocol as the plant growth regulators. From the docking outputs, compounds with the most favorable binding energy scores were selected within each ligand class, and the five strongest candidate interactions with microalgal proteins were identified. Predicted binding modes and interaction profiles were visualized and analyzed using BIOVIA Discovery Studio Visualizer [19].

2.4 Molecular Dynamics

Molecular dynamics (MD) simulations were performed using GROMACS 2024.1 [22] to evaluate the structural integrity and dynamic behavior of the protein–ligand complexes identified in docking (STGA/FabD, KASII/ZOG, FabG/STG, FATA/OROB, DHZMP/GPAT). For system preparation, hydrogen atoms were added to the protein structures, which were then solvated in a cubic water box with a 1.2 Å margin using the TIP3P water model. Na+ and Cl− ions were introduced to neutralize the overall charge. The AMBER99SB-ILDN force field was applied to proteins, while ligands were parameterized using the Generalized Amber Force Field 2 (GAFF2). Energy minimization was performed in two steps: 5000 iterations of the steepest descent algorithm followed by 5000 iterations of the conjugate gradient method, with a uniform 10 Å cutoff for van der Waals interactions. The equilibration phase consisted of a 1 ns run under constant volume and temperature (NVT) conditions at 298.15 K. Production simulations were then carried out for 100 ns under constant pressure and temperature (NPT) conditions, employing a 2 fs time step and periodic boundary conditions. Post-simulation trajectory analyses included root mean square deviation (RMSD), root mean square fluctuation (RMSF), hydrogen bonds profiles, radius of gyration (Rg), and solvent-accessible surface area (SASA). MATLAB [23] was used for visualization of the results.

2.5 Density Functional Theory (DFT)

Density functional theory (DFT) calculations were performed for the prioritized ligands using Gaussian 09 [24]. The geometries were optimized using the B3LYP functional and 6-311G (d, p) basis set [25]. Results were analysed with GaussView 5.0 [26]. Electronic descriptors were derived, including frontier molecular orbitals (FMO), molecular electrostatic potential (MEP), and global chemical reactivity parameters. These included total energy, dipole moment, HOMO and LUMO energies, energy gap (∆E), absolute hardness (η), global softness (σ), chemical potential (μ), electronegativity (χ), ionization energy (IE), electron affinity (EA), and electrophilicity index (ω) [27]. The minimum-energy conformer obtained from docking (PDB format) was used as the input structure for the DFT calculations. These descriptors provide supplementary electronic characterization of the ligands under gas-phase conditions and are intended to highlight candidate charge-transfer properties and kinetic stability.

3 Results and Discussion

3.1 Overview from Virtual Screening

Structure-based virtual screening was performed to prioritize candidate binders of microalgal enzymes involved in lipid metabolism. For each compound, the best docking mode was retained. Fifty-eight molecules showed predicted binding energies below −5.0 kcal/mol, consistent with candidate interactions at the active sites of the five target enzymes. These compounds span eight phytohormone families: strigolactones, cytokinins, brassinosteroids, gibberellins, auxins, jasmonates, abscisic acid, and salicylic acid. In contrast, betaines, polyamines, and ethephon exhibited weaker predicted affinities (above −5.0 kcal/mol) and are considered less promising (Table S3, Supplementary Data). This observation aligns with previous reports indicating that ethephon treatment did not alter lipid content or gene expression in microalgal biomass [28].

Twenty-two molecules were prioritized across phytohormone classes, including strigolactones (OROB, STG, 5-DS, STGA), cytokinins (TZROG, ZOG, TZ9G, DHZMP), brassinosteroids (CS, EPCS, BL, and EPBL), gibberellins (GA7, GA3, GA6, GA4), auxins (IBA, IPA2, PAA, OXIAA), as well as ABA and SA. These compounds represent the chemical diversity of phytohormones with predicted candidate interactions at microalgal lipid metabolism enzymes. This observation is consistent with physiological studies reporting that exogenous 24-epibrassinolide was associated with lipid accumulation, a higher proportion of C16:C18 fatty acids, and upregulation of lipid-related genes such as ACP thioesterase (FATA) in Botryococcus braunii B12 [29], biotin carboxylase (BC), and fatty acid desaturase (FAD) in Chlorella pyrenoidosa ZF [28], and ω-3 fatty acid desaturase (FAD) in Chlorella vulgaris strains [30]. Previous studies have reported that GA3 treatment in Botryococcus braunii B12 was associated with increased lipid accumulation and upregulation of genes related to fatty acid biosynthesis, including ACP, BC, FAD, KAS, and MCTK [29]. In Scenedesmus obliquus, treatments with BAP and GA under nitrogen limitation were reported to be associated with increased lipid accumulation and upregulation of genes involved in fatty acid and TAG biosynthesis, including ACP, SAD, FATA, and DGAT [31]. Furthermore, treatments with IAA and ABA in microalgal cells were reported to be associated with ~20–21% increases in lipid content, together with upregulation of fatty acid biosynthesis genes such as ACP, MCTK, and FATA [32]. Similarly, salicylic acid application in Phaeodactylum tricornutum was associated with an ~88% increase in lipid content [33].

3.2 Ligand Receptor Affinity

According to the docking results, STGA, ZOG, STG, OROB, and DHZMP emerged as the top-scoring ligands, with predicted binding energies ranging from −8.3 to −11.3 kcal/mol. These strigolactone and cytokinin derivatives exhibited stronger predicted affinities than several reference inhibitors and substrate analogues, consistent with candidate interactions at catalytically relevant regions of FabD, KASII, FabG, FATA, and GPAT (Table 2). The binding poses suggest possible competition with reference ligands, pending experimental validation. As these enzymes mediate key steps in fatty acid and triacylglycerol (TAG) biosynthesis, the prioritized ligands represent hypothesis-generating candidates for further study. This observation is consistent with previous reports in which strigolactones and cytokinins were associated with enhanced lipid accumulation and upregulation of genes involved in carbon fixation and lipid synthesis, including rbcL, accD, KAS III, DGAT, PEPC, ME [31,34,35], as well as ACP, SAD, FATA, FAD, MCTK, and DGAT [28,29,30,31]. These parallels support the prioritization of the identified compounds while underscoring the need for experimental validation to determine their biological relevance.

3.3 Ligand Receptor Interactions

The specific interactions between the selected ligands and microalgal enzymes were analyzed, with binding residues summarized in Table S4 (Supplementary Data). For FabD, the strigolactones, (+)-strigyl acetate (STGA) was predicted to to form three hydrogen bonds (ASN152, GLN156, and LYS182) and two alkyl contacts (LEU181 and LYS182) within the active site (Fig. 1A). These residues are located near conserved positions (ARG, SER, HIS, GLN, and LEU) that contribute to FabD catalytic function [14,36]. The predicted STGA binding mode partially overlaps with that of the natural substrate malonyl-CoA, which engages residues such as GLN10, GLN13, HIS87×2, ASN150, GLN156, SER12, GLY11, and SER88, and also shares contacts with the standard inhibitor corytuberine. Such overlapping interactions suggest possible competition with malonyl-CoA or reference inhibitors at the initiation site of FAS II, pending experimental validation.

In the predicted complex between cis-zeatin-O-glucoside (ZOG) and β-ketoacyl-acyl carrier protein synthase II (KASII), the ligand was predicted to form eight hydrogen bonds: five conventional interactions (ILE273, ASN407, ASN314, and ALA209×2) and three carbon-hydrogen contacts (THR274, PHE401, and ARG210) (Fig. 1B). These residues are situated in the immediate vicinity of the key catalytic quadruplet (CYS172, HIS307, LYS339, and HIS344). PHE401 has been described as a “gatekeeper” residue that undergoes conformational changes to regulate access to the binding pocket. The active site crevice, defined by a network of residues, including ASP269, THR274, PRO276, GLY280, ALA283, THR309, THR311, ALA313, ASN314, THR317, GLU318, GLY402, PHE403, GLY405, HIS406, and ASN407 [15], provides a candidate environment for ZOG interactions. Compared to with the substrate analogue PN7 and the inhibitor platensimycin, which engage catalytic and elongation residues more directly, ZOG appears to interact less directly with the elongation site. This binding pattern is consistent with a candidate interaction near residues involved in the first elongation step of type II fatty acid synthesis [37]. Experimental validation will be required to determine whether such interactions translate into activation, inhibition, or allosteric modulation.

In the predicted complex between 3-oxoacyl ACP reductase (FabG) and (+)-strigol (STG), hydrogen bonds were observed with ARG18, TYR158, LYS162, and GLY95, together with hydrophobic interactions involving ALA40 and ALA94 (Fig. 1C). These residues are located near the conserved catalytic triad (SER, TYR, LYS) and an asparagine residue that contribute to FabG activity [36,38,39]. The predicted STG binding profile partially overlaps with that of the cofactor NADPH and the standard inhibitor luteolin, which also engage residues at the cofactor-binding region. This overlap suggests a candidate interaction consistent with possible competition at the NADPH-binding site. Whether such binding events result in activation, inhibition, or allosteric modulation remains to be determined through enzymatic validation.

In the predicted complex between (+)-orobanchol (OROB) and acyl-ACP thioesterase (FATA), an α/β hydrolase involved in terminating fatty acid elongation, three hydrogen bonds were observed (ALA137, ARG182×2) together with hydrophobic contacts involving ALA121 and TRP202 (Fig. 1D). These predicted interactions resemble those reported for the substrate analogue oleic acid and the inhibitor WP2, which also engage residues such as TRP202, ALA121/125/137, PHE147, and ARG182. The overlap in hydrophobic contacts suggests that OROB may share candidate binding features with both substrate and inhibitor at the C16–C18 chain-termination site. This observation is consistent with the quantitative model proposed by Jing et al. [40], who identified 22 amino acid residues mainly located within the N-terminal hot-dog fold of the predicted CvFatB2a structure, including ARG124, ARG125, ALA63, and ALA192, which are likely involved in substrate specificity through modulation of enzyme–ACP interactions. Further experimental studies will be required to clarify whether such binding corresponds to functional modulation of FATA activity.

Finally, glycerol-3-phosphate acyltransferase (GPAT), which catalyzes the first rate-limiting step of de novo glycerolipid synthesis [41], was predicted to interact with DHZMP. The binding profile included an electrostatic interaction with ASP253, ten hydrogen bonds (SER193, LYS194, VAL191, GLU144, GLY235×3, ALA232, GLY170, and HIS141), and five hydrophobic interactions (ALA169, PRO147, VAL173, LEU188, and TYR167) (Fig. 1E). These predicted contacts overlap with catalytically relevant regions of GPAT, involving ASP253, LYS194, and GLU144. Previous studies [18] reported similar binding sites in MiGPAT1, where ACP interacts via hydrogen bonds and electrostatic forces to stabilize the acyl intermediate. The observed DHZMP binding pattern also partially resembles that of sn-glycerol-3-phosphate and the inhibitor FSG67, consistent with a candidate interaction at the sn-1 acylation site. Additional in vitro assays will be needed to evaluate if the predicted binding translates into functional modulation of GPAT activity.

Figure 1: Molecular docking of STGA/FabD (A), ZOG/KASII (B), STG/FabG (C), OROB/FATA (D), and DHZMP/GPAT (E) complexes and the interaction residues.

3.4 Molecular Dynamics

Molecular dynamic (MD) simulations were performed to explore the conformational stability and dynamic behavior of the protein–ligand complexes [42]. Stability was assessed using RMSD, RMSF, hydrogen bond profiles, radius of gyration (Rg), and solvent-accessible surface area (SASA). RMSD provides an overview of conformational stability, RMSF highlights residue-level flexibility, hydrogen bonds contribute to binding persistence, Rg reflects structural compactness, and SASA indicates solvent exposure [43]. For the STGA/FabD complex (Fig. 2), the backbone RMSD remained within ~0.1 nm, whereas the complex RMSD stabilized between 0.5–0.7 nm. The ligand RMSD was below 0.2 nm, and RMSF values were generally <0.3 nm. Stable Rg (1.37–1.67 nm) and SASA (125–140 nm2), together with 1–3 persistent hydrogen bonds, are consistent with a compact structure and limited conformational rearrangement.

Figure 2: MD simulation of STGA/FabD complex during 100 ns (RMSD, RMSF, H-bonds, Rg, and SASA profiles).

In the ZOG/KASII complex (Fig. 3), the backbone RMSD stabilized between 0.1–0.25 nm, with ligand RMSD <0.3 nm. The complex RMSD plateaued at 0.6–0.88 nm (20–70 ns) before converging to 0.4–0.6 nm. RMSF values were <0.3 nm, with localized flexibility in residues 100–150. Stable SASA (156–174 nm2), Rg (1.55–1.85 nm), and 1–11 hydrogen bonds suggest a well-maintained binding mode.

Figure 3: MD simulation of ZOG/KASII complex during 100 ns (RMSD, RMSF, H-bond, Rg, and SASA profiles).

For the STG/FabG system (Fig. 4), backbone RMSD remained within 0.1–0.3 nm, while ligand RMSD was <0.2 nm. A gradual increase in complex RMSD toward the end of the simulation suggested ligand repositioning. RMSF highlighted flexibility in loop regions (100–110 and 200–220). Variations in Rg (1.25–1.60 nm) and SASA (111–127 nm2) are consistent with partial ligand displacement, although 1–5 hydrogen bonds were maintained.

Figure 4: MD simulation of STG/FabG complex during 100 ns (RMSD, RMSF, H-bonds, Rg, and SASA profiles).

The OROB/FATA complex (Fig. 5) reached equilibrium rapidly, with protein and ligand RMSD stabilizing below 0.4 nm after 10 ns. RMSF analysis revealed expected flexibility in terminal and loop regions. Stable hydrogen bonding (1–4 interactions), SASA (~145–163 nm2), and Rg (1.4–1.8 nm) support a compact and equilibrated structure.

Figure 5: MD simulation of OROB/FATA complex during 100 ns (RMSD, RMSF, H-bonds, Rg, and SASA profiles).

Similarly, the DHZMP/GPAT complex (Fig. 6) exhibited stable dynamics, with RMSD values maintained between 0.15–0.2 nm over 100 ns. Ligand and protein RMSD stabilized at ~0.15 and ~0.3 nm, respectively. RMSF indicated limited flexibility confined to the N-terminal region, whereas active site residues remained stable. Persistent hydrogen bonding (1–10 interactions), along with stable Rg (1.5–2.0 nm) and SASA (172–190 nm2), confirmed structural integrity.

Figure 6: MD simulation of DHZMP/GPAT complex during 100 ns (RMSD, RMSF, H-bonds, Rg, and SASA profiles).

Overall, all complexes maintained stable binding modes with limited backbone deviations and persistent key interactions throughout the simulations. These findings support the structural plausibility of the docking poses; however, they reflect only conformational stability, and the potential effects on catalytic efficiency or pathway flux must be validated through enzymatic and cellular experiments.

3.5 Density Functional Theory

3.5.1 Frontier Molecular Orbital (FMO) Analysis

Frontier molecular orbital (FMO) analysis was performed to characterize the electronic properties of the candidate ligands. Total energy values (−1113.31 to −1806.39 Hartrees) reflect intrinsic electronic stability, with more negative values indicating greater stability. Dipole moments ranged from 4.08 to 9.07 Debye, suggesting variable polarity that may influence potential interactions. Higher HOMO energies are associated with enhanced electron-donating capacity, whereas lower LUMO energies indicate increased electrophilicity, polarization, and intramolecular charge transfer, which may contribute to candidate binding potential [43]. The HOMO–LUMO energy gap (2.96–5.22 eV) provides an index of reactivity and stability. Fig. 7 shows the atomic orbital components of these compounds. For example, ZOG exhibited a relatively large gap (5.22 eV), consistent with kinetic stability and a “hard” electronic character [44]. In contrast, DHZMP displayed a smaller gap (2.96 eV), consistent with a “soft” system that is more polarizable and reactive [45]. These descriptors provide a theoretical basis for prioritizing ligands as supplementary electronic characterizations, highlighting their ability to adapt electron density to enzyme active-site environments and potentially favor interactions with polar residues or intermediates. Their biological relevance remains qualitative and requires experimental validation.

Figure 7: Molecular boundary orbitals of STG, STGA, OROB, DHZMP, and ZOG.

3.5.2 Global Chemical Reactivity Descriptors (GCRD)

HOMO and LUMO energies are related to ionization potentials and electron affinities. The ionization energies (4.38–6.84 eV) indicate thermal stability, while the electron affinities (1.14–2.65 eV) suggest moderate electron-accepting capacity. Absolute hardness values (1.48–2.61 eV) reflect resistance to electron exchange, whereas global softness (0.19–0.33 eV−1) highlights the ease with which electrons can be redistributed under external perturbation [44]. The negative chemical potential (−2.90 to −4.73 eV) is consistent with electronic stability, while electronegativity values (2.90–4.73 eV) reflect the tendency to capture electrons for bond formation. Notably, strigolactones (STG, STGA, OROB) exhibited relatively high electrophilicity indices, suggesting a stronger tendency to accept electron density. This observation is consistent with candidate interactions at catalytically relevant regions, particularly nucleophilic residues [43] (Table 3). Experimental studies will be required to clarify whether such electronic properties correspond to functional modulation of enzyme activity.

3.5.3 Molecular Electrostatic Potential (MEP)

Molecular electrostatic potential (MEP) maps (Fig. 8) were generated to visualize electrophilic and nucleophilic regions predicted to participate in non-covalent interactions [27]. In strigolactones, negative potentials regions (red) localized on oxygen atoms of the cyclohexene and γ-butyrolactone rings suggest candidate hydrogen bonding or electrostatic contacts. For DHZMP, negative regions around the N–H groups of the purine ring (bicyclic A–B) and the O–P group of the side chain are consistent with predicted hydrogen bonding and electrostatic interactions with GPAT. These spatial distributions qualitatively support the interaction patterns observed in docking and MD simulations.

DFT-derived descriptors provide useful insights into the intrinsic electronic properties of ligands under gas-phase conditions and do not capture the full complexity of the biological environment. Therefore, these results should be interpreted as supplementary, rather than predictive of catalytic activity. Any mechanistic conclusions will require advanced modeling approaches, such as QM/MM calculations, together with experimental validation.

Figure 8: Molecular electrostatic potentials of STG, STGA, OROB, DHZMP, and ZOG.

4 Conclusion

This in silico study represents an exhaustive screening of sixty-five plant growth regulators against key microalgal enzymes (FabD, KASII, FabG, FATA, and GPAT) involved in lipid metabolism. Through a multi-layered computational approach, strigolactones and cytokinins‒specifically OROB, DHZMP, and STGA‒were prioritized as candidate binders for further investigation. Their predicted binding energies, dynamic stability profiles, and complementary electronic descriptors provide a computational rationale for possible interactions at catalytically relevant regions of microalgal lipid metabolism. These hypothesis-generating results require experimental validation through enzymatic assays, microalgal cultivation, lipid quantification, and gene expression analyses to determine their biological relevance. Such efforts will be critical in guiding the rational design of biostimulation strategies aimed at optimizing microalgal lipid productivity and strengthening the economic viability of sustainable biofuels.