Structure of Genes Encoding Oxidosqualene Cyclases – Key Enzymes of Triterpenoid Biosynthesis from Sea Cucumber Eupentacta fraudatrix

1. Introduction

The Far-Eastern sea cucumber Eupentacta (=Cucumaria) fraudatrix is a producer of promising biologically active substances; in particular, a series of triterpenoid glycosides (saponins, so-called as cucumariosides) [1,2], exhibiting ichthyotoxic, antitumor, hemolytic, antifungal and other pharmacological effects [3,4,5,6]. It has been found that in sea cucumber E. fraudatrix all free, resulted by de novo biosynthesis, and are mainly of parkeol type with Δ⁹⁽¹¹⁾ bond while saponins are predominantly formed from lanosta-7,24-dienol (lanostadienol) type with Δ⁷⁽⁸⁾ bond in the triterpene scaffolds [7]. Besides a membrane stabilizing function, unusual sterols of E. fraudatrix were found to reduce the membranolytic action of own triterpenoid glycosides [7,8].

To date, metabolic pathways of triterpenoid biosynthesis have been partially characterized in a limited number of sea cucumber (holothurian) species [9,10,11]. The biosynthetic precursor of triterpenes, sterols and triterpene glycosides is 2,3-oxidosqualene, which is cyclized by a membrane oxidosqualene cyclase (OSC; EC 5.4.99.7). Among representatives of echinoderms such as sea stars (starfish) and sea urchins, one type of OSC (named lanosterol synthase, LSS) is encountered, catalyzing the cyclization of 2,3-oxidosqualene into lanosterol. In sea cucumbers, two OSC genes have been reported in four species Apostichopus japonicus, Stichopus horrensis, Stichopus chloronotus and Parastichopus parvimensis (now Apostichopus parvimensis) [10,11,12]. Based on its expression in lanosterol-deficient yeasts and the conversion of 2,3-oxidosqualene into parkeol or lanostadienol, these genes were named as parkeol synthases (PSs) and lanostadienol synthases (LDSs), correspondingly [10]. Phylogenetic analysis revealed PSs and LDSs grouped together forming two distinct clades on the isolated branch, which stands alone from other animal OSCs. That suggests evolution of two divergent OSCs from an ancestral sea cucumber LSS by gene duplication and neofunctionalization [10].

The previously published data make an important contribution to understanding of triterpenoid biosynthesis in sea cucumbers, concern to only the family Stichopodidae (order Synallactida). Herein, we present the results of the transcript and gene structures determination of two OSCs from E. fraudatrix, whose deduced amino acid sequences showed notable difference when compared with OSCs from Stichopodidae. The molecular modeling and docking analyses revealed that OSC1 might be parkeol synthase and OSC2 might be lanostadienol synthase. These findings, firstly acquired for family Sclerodactylidae (order Dendrochirotida), will bring some additional clarity on the origin and biosynthesis of triterpenoid glycosides in sea cucumbers.

2. Materials and Methods

2.1. Sample Collection and Nucleic Acids Extraction

The samples of sea cucumber E. fraudatrix were collected in Troitsa and Sobol Bays (Peter the Great Gulf, the Sea of Japan) during 2019 – 2024. Species identification was carried out by Dr. I.Y. Dolmatov (A.V. Zhirmunsky National Scientific Center of Marine Biology, FEBRAS, Vladivostok, Russia). Genomic DNAs and total RNAs were extracted from frozen or fresh E. fraudatrix tissues (body wall, muscle, gut) using TRIzol reagent (Invitrogen, USA) according to the manufacturer`s manual, except that the final RNA elution step was performed with HiDi Formamide (Thermo Fisher Scientific, Waltham, MA, USA) and stored in −70 ℃. Nucleic acids concentration and purity were determined on a NanoPhotometer P330 (Implen, Munich, Germany), and their integrity was assessed using electrophoresis in 1.2% agarose gel.

The study was carried out in accordance with the recommendations of the Convention on Biological Diversity and was approved by the Ethics Committee of the G.B. Elyakov Pacific Institute of Bioorganic Chemistry, FEBRAS, Vladivostok, Russia (Protocol No. 0037. 12 March 2021).

2.2. cDNA Synthesis, RACE and PCR Amplification

The synthesis of full-length-enriched double stranded cDNAs was carried out using total RNAs, pretreated with DNase I (Thermo Fisher Scientific, Waltham, MA, USA) and MINT Universal reagent kit (Evrogen, Moscow, Russia) in accordance with the manufacturer`s protocols. The 3`- and 5`- flanking sequences of OSC cDNAs were obtained by the Rapid Amplification of cDNA Ends (RACE) technique using the RACE primer set (Evrogen, Moscow, Russia) and gene-specific primers (Table 1). The primers were designed based on partial OSC sequences from E. fraudatrix transcriptome GHCL00000000 [13]. The RACE amplifications were conducted with Encyclo® DNA Polymerase (Evrogen, Moscow, Russia) according to the scheme (Figure S1). RACE-PCR-fragments were cloned and sequenced using ABI 3130xl Genetic Analyzer (Applied Biosystems, Hitachi, Tokyo, Japan). The cDNA sequences were verified by PCR amplification with OSC1_start/ OSC1_stop and OSC2_start/ OSC2_stop primers (Table 1) using Q5® High-Fidelity DNA Polymerase (New England Biolabs, Ipswich, MA, USA) under following temperature conditions: 98°C for 30 s, 35 cycles of 98°C for 10 s, 62°C for 15 s, and 72 °C for 70 s, and finally 72 °C for 300 s. Obtained full-length CDSs were cloned using the InsTAclone PCR Cloning Kit (Thermo Fisher Scientific, Waltham, MA, USA) and positive clones were sequenced with gene-specific primers using SeqStudio™ Genetic Analyzer (Thermo Fisher Scientific, Waltham, MA, USA). The obtained sequences were deposited into GenBank databases under the accession numbers OR711403.1 and OR725688.1 respectively for OSC1 and OSC2.

2.3. Gene Structure Determination

Amplification of OSC1 and OSC2 gene fragments from genomic DNAs with size be apart in 5 – 12 Kbp were performed with Long PCR Enzyme Mix (Thermo Fisher, Waltham, MA, USA) and gene-specific primers (Table 1) under following temperature conditions: 94℃ 1 min, 30 cycles of 94℃ 20 sec, 58℃ 30 sec, 68℃ 600 sec. Fragments with size larger than 12 Kbp were obtained with KOD One Master Mix (Toyobo, Osaka, Japan) according to the manufacturer`s protocol. The ONT DNA libraries were prepared using SQK-NBD114.24 (Oxford Nanopore Technologies, Oxford, UK) and sequencing with MinION Mk1b on flow cell FLO-MIN114 (Oxford Nanopore Technologies, Oxford, UK). The ONT raw reads were basecalled with Dorado v.0.7.1 (Oxford Nanopore Technologies, Oxford, UK). The BRAKER version 3.0.6 [14] was used to predict de novo OSC genes and untranslated regions (UTRs) using the both RACE-PCR and transcriptomic data of E. fraudatrix [13], protein dataset of Metazoa from OrthoDB version 11 [15], and proteins of Apostichopus japonicus [16], Holothuria leucospilota [17], Acanthaster planci (GCF_001949145.1) [18], Asterias rubens (GCF_902459465.1) [19], Strongylocentrotus purpuratus (GCF_000002235.5) [20], Lytechinus variegatus [21] and Patiria miniata (GCF_015706575.1) uploaded from NCBI [22] and Echinobase [23]. Protein containing sequences were extracted and processed through QIAGEN CLC Genomic Workbench, version 24.0.1 (QIAGEN, Aarhus, Denmark).

2.4. Sequence Analysis, Homology Modeling and Docking

Multiple sequence alignments of OSCs were performed using ClustalW implemented in MEGA 11, version 11.0.13 [24] or QIAGEN CLC Genomic Workbench, version 24.0.1 (QIAGEN, Aarhus, Denmark) and manually corrected. A maximum-likelihood (ML) phylogeny with 500 non-parametric bootstrap replicates was obtained under LG+G substitution model calculated in MEGA 11, version 11.0.13 [24].

Theoretical models of EfOSC1, EfOSC2, AjLAS1, and AjLAS2 were obtained using the MOE package, version 2018.01 (Molecular Operating Environment; Chemical Computing Group ULC, 1010 Sherbrooke St. West, Suite #910, Montreal, QC, Canada, H3A 2R7, 2018). For molecular docking, the models of OSCs complexes with parkeol or lanostadienol were optimized by forcefield Amber10:EHT, using the crystal structure of human OSC in complex with lanosterol (PDB code 1W6K) [25] as a template after removing water molecules. The MOE docking parameters were set using default values. Contact analysis and visualization of the results were carried out using Ligand interaction and Dock modules in the MOE 2018.01 program and Discovery Studio Visualizer v24.1.0.23298 (Dassault Systèmes: San Diego, CA, USA, 2024). Structure alignments were performed using QIAGEN CLC Genomic Workbench, version 24.0.1 (QIAGEN, Aarhus, Denmark).

3. Results and Discussion

3.1. E. fraudatrix Genome Codes Two Different OSCs

3.1.1. cDNA Determination of OSC mRNAs

To design gene-specific OSC primers for RACE reactions, the draft transcriptome of E. fraudatrix [13] was used for a BLASTX search with the human LSS sequence as a query. The strategy for determination of OSC encoding transcripts is presented in Figure S1. As a result of a series of 5`- and 3`-RACE reactions, different RACE-products were obtained depending on tissues E. fraudatrix cDNA libraries; the 5`-fragments were about 600 bp (body wall) and 400 bp (gut) and the 3`-fragments were about 600 bp (body wall) and 1000 bp (gut). It means that at least two OSC isoforms (OSC1 and OSC2) could be transcribed from the sea cucumber E. fraudatrix genome. The RACE sequences were further used to design primers flanking OSC1 CDS (pair primers OSC1_start/ OSC1_stop) and OSC2 CDS (pair primers OSC2_start/ OSC2_stop). After amplification, cloning and sequencing, the OSC1 and OSC2 CDSs derived from the body wall and gut cDNA libraries had the same size of 2163 bp.

Totally, after assemblies of CDS and RACE fragments, the full-length OSC1 and OSC2 cDNAs (= open reading frame, ORF) besides the CDS include 5`-untranslated regions (UTRs) of at least 105 bp and 299 bp, and 3`-UTRs of at least 141 bp and 801 bp, respectively. The OSC1 and OSC2 mRNAs encode a 721-amino-acid protein with calculated molecular weights of 82.531 and 82.502 kDa, respectively. Importantly, that the predicted 5’UTRs of OSCs genes derived from draft E. fraudatrix genome (unpublished data) were validated with presented RACE-PCR results where both 5`-UTRs started with G at position +1 and the first initial AUG codons with a different Kozak consensus sequence of -9GTGAGAGCG-1 for OSC1 mRNA and -9CTGCAAGGA-1 for OSC2 mRNA. Interestingly, the 5′-ends of mRNAs contained multiple ATG codons, however all but one translated in truncated proteins. However, the obtained 3'-UTRs sequences did not confirm the BREAKER prediction results, since Sanger sequencing reactions were interrupted by repeats and poly(T) regions. In the obtained OSC1 3`-UTR, a polyadenylation signal site was not detected, whereas in the OSC2 3`-UTR it was located at +3004AATAAA+3009. Comparative cDNA analysis showed that the 5`-UTRs and 3`-UTRs did not share any significant similarity, however OSC1 and OSC2-encoding sequences shared of 79.96% similarity.

3.1.2. Sequence Analyses of OSC1 and OSC2 Proteins

Based on Interpro search results, the deduced amino acid sequences of E. fraudatrix OSC1 and OSC2 belong to terpene synthase (IPR002365, https://www.ebi.ac.uk), which combines the evolutionary related lanosterol synthase (EC 5.4.99.7), cycloartenol synthase (EC 5.4.99.8) and hopene synthase (EC 5.4.99.) [26]. In addition, BLASTP/CDD searches showed that the OSCs belong to ISOPREN_C2_like superfamily, Class II terpene cyclases, including squalene cyclase and 2,3-oxidosqualene cyclase [27]. These enzymes are integral membrane proteins, which are responsible for converting linear triterpenes such as squalene or oxidosqualene into cyclic triterpenes such as hopene, lanosterol or others.

Comparison of the primary structures of OSCs and LSSs from echinoderms (Table 2) showed that the E. fraudatrix OSCs sequences shared relatively high identity (65–71%) and similarity (79–83%) with other holothurian OSCs. They demonstrated less than 58% identity and 73% similarity to LSSs from the sea urchin and starfish that might indicate evolutionary divergence in their sterol biosynthesis pathways. The primary structures of holothurian OSCs had much lower identity (53–57%) and similarity (69–73%) to human LSS.

Comparison of the E. fraudatrix OSCs sequences with those of OSCs from sea cucumbers, sea urchin, starfish and human has carried out by the multiple alignment of amino acid sequences (Figure 1). Like other OSCs, E. fraudatrix OSC1 and OSC2 contain the major functional domains (Figure 1): four or five QW motifs to stabilize the enzyme structure [28,29], highly conserved DTTAE domain to bind and protonate a substrate [30,31,32,33], and LWIHCR domain to provide the cyclization of a substrate [34,35]. Interestingly, the DTTAE motif is conserved in both E. fraudatrix OSCs, whereas in Apostichopus and Stichopus AjLAS1, ApPS and ShOSC1 have DTTAE, and AjLAS2, ApLDS and ShOSC2 have DTSAE. Despite the established difference in catalytic function of AjLAS1, ApPS and AjLAS2, ApLDS [10,12], the substituted amino acids share similar properties, and these enzymes use the same substrate oxidosqualene. Moreover, LSSs presented on the Figure 1 have DCTAE motif.

Notably, the sequence identity and even the sequence similarity alone did not help to identify the OSC orthologues (Table); E. fraudatrix OSC1 and OSC2 had higher sequence identity to each other (77.5%) than to A. japonicus LAS1 or LAS2 (67.3–70.6%) (Table 2). Because in our case the sequence identity between holothurian paralogues (E. fraudatrix, family Sclerodactylidae) is greater than between holothurian orthologues (families Sclerodactylidae and Stichopodidae), it is very likely that the similarity of specific amino acid residues between the active sites should be more useful in identifying orthologues which perform a similar function by a definition. The active site of OSCs is mainly represented by highly conserved amino acid residues. However, superimposition of the HsLSS, ApLDS and ApPS homology models performed by Thimmappa et al. [10] revealed that the amino acid residue at position 444 alone can distinguish PSs (436L) and LDSs (444Q) from LSSs (444F). Comparison of the active sites showed that EfOSCs do differed at this residue, which resulting in the grouping of EfOSC1 with AjLAS1 (=ApPS) at 436L and EfOSC2 with AjLAS2 (=ApLDS) at 439Q (Figure 1).

3.2. Confirmation of Key L and Q Residues Distinguishing Parkeol and Lanostadienol Synthases by Molecular Docking

A comparative modeling approach was taken to infer the function and orthology of E. fraudatrix OSC1 and OSC2 from their structures and determine the specific amino acid residues required to form Δ9(11) or Δ7(8) bonds. Importantly, two homologues AjLAS1 and AjLAS2 recently identified in the A. japonicus genome, and characterized by yeast heterologous expression showed that AjLAS1 catalyzes the formation of Δ9(11) bond (parkeol), whereas AjLAS2 synthesizes lanostadienol with Δ7(8) bond [11].

In addition, corresponding mutations were later introduced in A. japonicus PS (=AjLAS1) as well as in sea urchin S. purpuratus LSS, which verified the key role of leucine (436L) in the biosynthesis of parkeol and glutamine (444Q) in the biosynthesis of lanostadienol [10]. Namely, Gil77 yeast cells with A. japonicus PS^L436F and PS^L436Q mutants synthesized lanosterol and lanostadienol, respectively, in addition to parkeol. The Gil77 yeast cells with the mutant of S. purpuratus LSS^F440L showed detectable levels of parkeol, whereas Gil77 cells with the mutant of S. purpuratus LSS^F440Q synthesized detectable levels of lanostadienol.

To compare the putative active sites of EfOSC1 and EfOSC2 to those of the functionally characterized AjLAS1 (parkeol) and AjLAS2 (lanostadienol), their structural models were generated by homology modeling using human LSS [PDB 1W6K] as a prototype [25]. All models showed root-mean-square deviations of Сα-atom of less than 1Å (0.18Å), indicating very similar overall structures with the prototype. The superposition of the obtained OSCs structures revealed TM-scores of 0.977 to 0.986 and RMSD values over all atoms of 0.928 to 1.392 Å (Table S1, Figures S2, S3). Modeling of enzyme-product docking with MOE S-score lower than −11.5 identified pairwise differences in the binding sites; docking of parkeol into the EfOSC1 or AjLAS1 as well as lanostadienol into the EfOSC2 or AjLAS2 resulted in the local structural similarities despite the lack of overall sequence identities. The 2D structures of EfOSCs and AjLASs with active site residues within 5Å around the products are shown in Figure S4. In contrast, docking of lanostadienol into the EfOSC1 or AjLAS1 and parkeol into the EfOSC2 or AjLAS2 led to similar patterns of unfavorable H-bonds (Figure S5).

Based on the resulting 3D models (Figure 2), the location of residues of leucine (436L/435L) in EfOSC1/AjLAS1 and glutamine (439Q/444Q) in EfOSC2/AjLAS2 as well as threonine (448T/447T or 451T/456T), tyrosine (Y496/Y495 or Y498/Y503) and phenylalanine (F688/F690 or F691/F696) ones near the B ring of lanostadienol (Figure 2c,d) or C ring of parkeol (Figure 2a,b) is characterized by high structural similarity, suggesting that these particular residues 436L and 439Q may determine product specificity in the EfOSCs.

Thus, based on the results of molecular docking, it can be inferred that in sea cucumber E. fraudatrix OSC1 should cyclize 2,3-oxidosqualene into parkeol and be an orthologue of parkeol synthase AjLAS1 (=AjPS), while OSC2 should cyclize it into lanostadienol and be an orthologue of lanostadienol synthase AjLAS2 (=AjLDS). Further biochemical studies of the OSCs proteins and their mutants need to be carried out to prove the putative function of these enzymes.

3.3. E. fraudatrix OSC1 and OSC2 form a Phylogenetically Distinct Branch

To evaluate the evolutionary relationships of E. fraudatrix OSC1 and OSC2 in Echinodermata, phylogenetic trees were constructed using the ML method based on aligned amino acid sequences of known OSCs (PSs, LDSs, and LSSs) from Metazoa (Figure S6). Using OSCs from Dictyostelium discoideum (Amoebozoa, sister clade to Obazoa) and Saccharomyces cerevisiae (Holomycota, sister clade to Holozoa) as outgroups [36], OSCs from Metazoa formed two super-clades (Figure 3). The first clade exclusively represented holothurian OSCs (E. fraudatrix, A. japonicus, A. parvimensis and S. horrens), and the second clade included LSSs from Porifera, Placozoa, and Bilateria as well as other echinoderms such as starfish and sea urchins. In the starfish and sea urchin genomes, a single OSC gene was identified, which, like in all other animals, encodes lanosterol synthase with F444. Based on our phylogenetic analysis, holothurian OSCs were grouped into clades according to their taxonomic positions; the first clade consisted of two paralogues OSC1 and OSC2 from E. fraudatrix (family Sclerodactylidae) and the second clade consisted of two subclades corresponding PSs and LDSs from members of the genera Apostichopus and Stichopus (family Stichopodidae) (Figure 3a). In case of exclusion of the holothurian OSCs from the alignment used for tree construction, the resulting tree reflected the taxonomic relationships inferred from their genome sequences (Figure 3b).

It is noteworthy that on the genomic tree [21], holothurians and sea urchins were grouped more closely, while sea stars are situated on the outer branch of Echinodermata. In contrast, in the OSC tree constructed, sea urchins clustered with sea stars, whereas holothurians formed an outer branch for the Metazoa clade with high branch node supporting. This suggests that despite the phylogenetic proximity of holothurians to other echinoderms, they demonstrate a distinct independent evolutionary lineage of OSCs driven by duplication events that provide the biosynthesis of unusual sterols and triterpenoid glycosides.

Together, despite that E. fraudatrix OSCs formed a distinct branch rather than complementing orthologous groups of LDSs and PSs, we consider that the key residue in the position 444 (436L or 439Q) may determine not only the specific cyclization of the substrate, but also an attribution of the sequences to a certain orthologous group. The genetic distances observed between E. fraudatrix OSCs are significantly large, suggesting that their evolutionary divergence was driven by an acquired end-product specificity (in contrast to A. japonicus, E. fraudatrix synthesizes parkeol sterols and lanostadienol glycosides).

3.4. Gene Structure Determination and Analysis of OSC Genes

The EfOSC1 and EfOSC2 encoding gene sequences were obtained by a series of PCRs using gene-specific primers designed based on the cDNA sequences (Table 1). As a result of PCR optimization, fragments ranging from 4 to 18 kb were obtained. The resulting amplicons were sequenced using Oxford Nanopore Technologies, assembled using CLC Genomic Workbench 24.0.1 and manually corrected. The final OSCs gene sequences had a length of about 30-40 kb for OSC1 and 22-25 kb for OSC2. In addition, OSC1 and OSC2 cDNA sequences were mapped on E. fraudatrix draft genome (unpublished data) using BLAST, and regions with best hits were extracted. A total of four sequences were taken for each OSC gene (two from the amplicons assembly and two from the draft genome sequence), which were used for comparison with cDNA sequences and further haplotype analysis. After alignment with cDNA sequences, it became clear that the range of gene lengths was due to significant differences in intron sizes. The exon-intron structures of E. fraudatrix OSCs genes were comprised of 19 exons and 18 introns (Figure 4a) with conservation of the GT/AG splice sites in both genes. The OSC1 5'UTR was located in exon 1, while the OSC2 5'UTR occupied exon 1 and partially exon 2. The last codons of CDSs and 3’ UTRs were aligned on the same 19th exon for both genes. The lengths of the corresponding protein-coding exons of OSC1 and OSC2 were equal, whereas the lengths of the corresponding introns varied significantly and ranged from 0.4 to 7.9 kb. Moreover, the lengths of some introns varied among haplotypes; for example, the size of OSC1 intron 2 ranged from 1.1 to 4.4 kb (Figure 4a).

Although the OSC1 and OSC2 transcript sequences show a high degree of similarity (Figure 4b, 1), alignment of the gene sequences using BLAST Global Alignment revealed an almost complete absence of aligned regions (Figure 4b, 3), indicating an exceptionally high degree of intronic divergence between E. fraudatrix OSC1 and OSC2 genes. This suggests a significant evolutionary rate for intronic regions while preserving overall sequence similarity at the transcript and protein levels (Figure 4b, 1,2). When OSC1 and OSC2 haplotypes were aligned (Figure 4c, d), a significant degree of polymorphism was found in OSC1 introns, whereas OSC2 introns showed greater conservation. The high polymorphism and larger size of the OSC1 gene suggest that this gene may be an ancestor of the OSC2 gene.

4. Conclusions

Our studies indicate that sea cucumber E. fraudatrix (family Sclerodactylidae, order Dendrochirotida) has two OSC genes, as well as representatives of the genera Apostichopus and Stichopus (family Stichopodidae). Based on the results of sequence analysis and molecular docking, we suggest that OSC1 with residue L436 seems like to be parkeol synthase, while OSC2 with residue Q439 appears to be lanostadienol synthase, and the above allows us to assume that these residues are key in identifying orthologous groups. We believe that if E. fraudatrix OSCs form their own separate phylogenetic clade, then this evolutionary divergence between holothurians of families Sclerodactylidae and Stichopodidae was probably driven by difference in the use of triterpene scaffolds for sterol and glycoside biosynthesis. Considering the transcript and gene structures of E. fraudatrix OSCs haplotypes, we suggest that OSC1 gene can be an ancestor of the OSC2 gene. Taken together, these findings, firstly acquired for family Sclerodactylidae (order Dendrochirotida), bring some additional clarity on the origin and biosynthesis of triterpenoid glycosides in sea cucumbers.