Genomic Insights into Fertilization: Tracing PLCZ1 Orthologs Across Amphibian Lineages

Abstract

Fertilization triggers a cascade of events, including a rise in egg cytosolic calcium that marks the onset of embryonic development. In mammals and birds, this critical process is mediated by the sperm-derived phospholipase C zeta (PLCζ), which is pivotal in releasing calcium from the endoplasmic reticulum in the egg and initiating embryonic activation. Intriguingly, Xenopus laevis, a key model organism in reproductive biology, lacks an annotated PLCZ1 gene, prompting questions about its calcium release mechanism during fertilization. Using bioinformatics and RNA sequencing of adult X. laevis testes, we investigated the presence of a PLCZ1 ortholog in amphibians. While we identified PLCZ1 homologs in 25 amphibian species, including 14 previously uncharacterized orthologs, we found none in X. laevis or its close relative, Xenopus tropicalis. Additionally, we found no compensatory expression of other PLC isoforms in these species. Synteny analysis revealed a PLCZ1 deletion in species within the Pipidae family and another intriguing deletion of potential sperm factor PLCD4 in the mountain slow frog, Nanorana parkeri. Our findings indicate that the calcium release mechanism in frog eggs involves a signaling pathway distinct from the PLCζ-mediated process observed in mammals.

Introduction

For most animals, fertilization triggers a surge in the calcium levels within the egg cytoplasm that begins near the site of sperm entry and permeates the entire egg (Steinhardt et al. 1977; Gilkey et al. 1978). This calcium catalyzes essential processes that prevent polyspermy and initiate embryonic development (Denninger et al. 2014; Swann and Lai 2016). Among the hypothesized mechanisms for this signal in mammals (Saunders et al. 2002; Nomikos et al. 2017) and birds (Coward et al. 2005) is the release of the sperm-derived enzyme phospholipase C zeta (PLCζ) into the egg. PLCζ cleaves the phospholipid phosphatidylinositol 4,5-bisphosphate (PIP₂) to release soluble inositol trisphosphate (IP₃) within the egg, which then initiates the release of calcium from the endoplasmic reticulum (ER) (Cox et al. 2002; Saunders et al. 2002; Yoneda et al. 2006; Mizushima et al. 2007; Ito et al. 2008; Ross et al. 2008; Bedford-Guaus et al. 2011; Mizushima et al. 2014). While ER-derived calcium is pivotal for egg activation in most animals (Crossley et al. 1988; Whitaker 2006), it remains uncertain whether the PLCζ-anchored mechanism is similarly widespread.

The African clawed frog (Xenopus laevis) is a long-standing and fundamental model organism for fertilization research. Similar to mammals and birds, fertilization in X. laevis eggs induces an increase in cytosolic calcium (Grey et al. 1982; Kline 1988). This calcium surge triggers the fast and slow blocks to polyspermy (Wozniak, Phelps, et al. 2018; Wozniak and Carlson 2020) and sets embryonic development in motion (Whitaker 2006). Despite these similarities, no PLCζ-encoding gene (PLCZ1) has been identified in the genomes of X. laevis or the related Xenopus tropicalis. While PLC activation is integral for this process (Wozniak, Tembo, et al. 2018), our previous studies found no evidence of activation in the PLC isoforms (Plcg1, Plcb1, and Plcb3) expressed in X. laevis eggs during fertilization (Komondor et al. 2023). One possibility is that X. laevis sperm donates an unidentified soluble sperm factor like PLCζ to the egg, as occurs in mammals and birds. Notably, X. laevis sperm extracts introduced into mouse eggs produce a calcium rise comparable to that caused by chicken or cow sperm extracts (Dong et al. 2000). Using a comparative genomics and transcriptomics approach, we set out to identify potential sperm factors in Xenopus and other amphibian species with readily available genome or transcriptome assemblies.

Results

At the outset of our study, we searched for genes encoding PLCζ in amphibians. Our initial analysis revealed only 11 annotated amphibian PLCZ1 orthologs, making it the least represented among all tetrapod classes (Fig. 1a). To broaden our investigation, we employed tblastn using the protein sequence of PLCZ1 from the plains spadefoot toad (Spea bombifrons) as a query against the NCBI transcriptome shotgun assembly (TSA) database and other transcriptome databases (Table 1). Additionally, we searched for PLCZ1 in amphibian whole genome shotgun (WGS) assemblies, using miniprot (Li 2023) to identify open reading frames. We required a query coverage of greater than 30% to consider identified sequences for further analysis (Table 2). These comprehensive searches identified not only PLCZ1 but also PLCD1-4 orthologs.

We differentiated these orthologs through multiple sequence alignment (MSA) and phylogenetic analysis alongside all annotated amphibian PLC isozymes. When grouped by clade, these sequences appeared across all 3 amphibian orders: tailless frogs and toads (Anura), newts and salamanders (Urodela), and limbless caecilians (Gymnophiona). However, we were unable to recover PLCZ1 orthologs from many TSA datasets, likely because PLCZ1 expression is typically restricted to the testes (Cox et al. 2002; Nomikos et al. 2005). Despite the presence of testis tissue in several Xenopus transcriptome assemblies, no evidence of PLCZ1 was detected (Fig. 1a).

Our initial search for PLCZ1 uncovered many PLCD1-4 transcripts, which became apparent when we aligned the sequences to annotated amphibian PLCs. We similarly found PLCD1-4 genes in our search within WGS assemblies. Phylogenetic analysis distinguished PLCZ1 orthologs from PLCD1-4 with robust bootstrap support (Fig. 1b; supplementary fig. S1, Supplementary Material online), reinforcing our original classifications. However, there remained a concern regarding the potential misidentification of PLCD1-4 as PLCZ1 due to their pronounced similarities. Additionally, some PLCD1-4 transcripts might lack the sequence encoding the determinative pleckstrin homology domain due to RNA-sequencing errors (Wang et al. 2009).

To further validate our findings, we performed a separate search for PLCD1-4 in the transcriptome assemblies where we had previously detected PLCZ1. PLCD1 and PLCD3 were found in all queried species except the Yenyuan stream salamander (Batrachuperus yenyuanensis), possibly explained by the low completeness (40.1%) of its assembly (Xiong et al. 2019). More surprisingly, we discovered species lacking a gene encoding PLCD4. This enzyme is expressed in multiple X. laevis tissues, including the testis (Bowes et al. 2010; Session et al. 2016), and is integral to the acrosome reaction in mammalian sperm (Fukami et al. 2001; Table 1).

To understand how anurans lacking PLCZ1 and PLCD4 diverged from related species, we examined their conserved synteny, which refers to the conservation of gene blocks on chromosomes across different species. We identified a shared synteny block containing PLCZ1 in the common frog (Rana temporaria), common toad (Bufo bufo), and the plains spadefoot toad (S. bombifrons). Using gene anchors from this block, we compared these species and frogs from the Pipidae family, including 3 Xenopus species and the Congo dwarf clawed frog (Hymenochirus boettgeri) (Fig. 2a). For assemblies lacking annotations, we used miniprot (Li 2023) to identify genes (Table 3). Our data showed that every available assembly from Pipidae lacks PLCZ1, indicating a gene deletion in that family's ancestral lineage.

We also reanalyzed species from our ortholog search to locate a syntenic block for PLCD4. In the mountain slow frog (Nanorana parkeri), this block revealed a gene deletion, although it does not occur in other members of its family (Hime et al. 2021; Fig. 2b).

We performed a transcriptomics analysis on X. laevis to assess the expression of various potential sperm factors and to detect potential unannotated PLCZ1 orthologs. Our best understanding of PLCζ function comes from studies in mammals, where it is exclusively found in sperm (Session et al. 2016). Mature sperm are transcriptionally quiescent (Kierszenbaum and Tres 1975), rendering them unsuitable for transcriptomics. Therefore, we opted to perform RNA-seq on whole testes. While we acknowledge that not all RNA in the testes translates to proteins in sperm, RNA encoding sperm proteins should be present within the testis. For our analysis, we extracted RNA from the testes of 3 adult X. laevis males. Following library preparation and sequencing, we produced an average of 66,723,678 contiguous RNA-seq reads per biological replicate, with 48% to 51% aligned to the X. laevis genome. We observed at least 1 transcript in 1 testis from an individual frog for 33,377 genes. We performed gene ontology (GO) analysis on 6,000 of the most highly expressed genes among these testis-derived transcripts. Of these, 48.3% were identified with GO biological process terms, many of which were essential for general cell maintenance (e.g. translation, protein folding, and cell cycle) (Table 4). Some were also involved in specialized processes such as spermatogenesis. Our overall results were highly correlated (supplementary fig. S2, Supplementary Material online) and matched a separate dataset (Fig. 3a). In our dataset, we identified several essential sperm transcripts, including deleted in azoospermia-like (dazl) (Houston and King 2000) and the testes-specific histone (h1-10) (Shechter et al. 2009; Oikawa et al. 2020). Other enriched transcripts, such as astl2b, rflcii, and ribc1, are known to be specifically expressed in X. laevis testis (Session et al. 2016) or as proteins in sperm (Sakaue et al. 2010; Ferlin et al. 2012; Agarwal et al. 2016; Maccarinelli et al. 2017). Some of these, including rflcii, cfh, and tubb4b, are essential for fertility in other species (Table 4; Sakaue et al. 2010; Ferlin et al. 2012; Agarwal et al. 2016; Maccarinelli et al. 2017).

We searched the RNA-seq dataset for transcripts encoding sperm factors, beginning with PLC isoforms. We identified transcripts for PLCβ, PLCδ, and PLCγ (Fig. 3b). Of these, the plcd4 transcript was the most abundant in our dataset. Notably, plcd4 was also the most abundant PLC transcript in testis datasets from the inbred X. laevis J strain (Session et al. 2016) and the X. tropicalis Nigerian strain, neither of which possessed a PLCZ1 ortholog. By comparison, testes of the Oriental fire-bellied toad Bombina orientalis contained PLCZ1 transcripts but also substantially reduced levels of plcd4 compared with Xenopus.

In mice, we observed that Plcz1 was the predominant PLC isoform expressed in a dataset from 3-month-old mouse testes. However, Plcd4 levels were comparable to those found in Xenopus species (Fig. 3b; Huang et al. 2021). We also assessed the presence of wbp2nl (PAWP), a sperm head protein proposed to be involved in egg activation in both frogs (Aarabi et al. 2010) and mice (Wu et al. 2007). While mouse and Xenopus exhibited comparable wbp2nl transcript levels, B. orientalis lacked this transcript entirely (Fig. 3b).

To further examine the absence of PLCZ1 in Xenopus species, we evaluated the unmapped RNA from our dataset for potential PLCZ1 orthologs. We conducted a blastx search using the unmatched reads against a database of PLCZ1 sequences from mammals and amphibians. This search yielded 360 reads with an E-value <10⁻⁶. These reads were then aligned to a PLCZ1 protein sequence from the neotropical leaf frog, Phyllomedusa bahiana, which was not included in our initial blastx database. The alignment covered residues 170 to 654 of this ortholog. However, further evaluation showed the reads were a closer match to PLCD1-4 (Supplementary material online). This suggests the unmapped reads in our dataset likely originate from incomplete contigs of annotated genes.

Discussion

As in other species, fertilization in X. laevis results in a surge of cytosolic calcium in the egg, a process facilitated by PLC enzymes that produce IP₃ to induce calcium release from the ER (Runft et al. 1999; Runft et al. 2002). Various mechanisms have been proposed for initiating this cascade of events. In mammals and birds, it is believed that sperm deliver soluble PLCζ to the egg during fertilization, serving as the trigger that initiates the calcium surge in the egg cytoplasm (Swann and Lai 2016).

Because of its role in mammalian and avian fertilization, we were surprised to find no annotated ortholog of PLCZ1 in X. laevis. This prompted us to investigate, which other amphibians possess this enzyme. Using annotated genes and TSA, we successfully identified PLCZ1 orthologs in 25 distinct amphibian species (Fig. 1b). These samples provide broad taxon coverage for amphibians, especially considering that PLCZ1 transcripts are hypothetically constrained to the testes, a tissue type frequently absent from transcriptomics datasets.

Our inability to detect a PLCZ1 ortholog in either Xenopus species implies a possible gene deletion in the Xenopus ancestral lineage, as corroborated by our local synteny analysis (Fig. 2a). When we included H. boettgeri in the analysis, a PLCZ1 deletion in a Pipidae common ancestor became apparent. We attempted to establish this deletion's origin for all sequenced family members of Pipidae. However, we were unable to extend this analysis to Surinam toads (Pipa) due to incompleteness of the genomic assemblies available for these species.

We also noted the absence of PLCD4 in N. parkeri (Fig. 1d). Local synteny comparison confirmed a PLCD4 deletion in this species, which appears to be a recent event since another Ranidae species, R. temporaria, maintains a functional version of this gene (Ma et al. 2018). We observed 3 additional species missing PLCD4. The Yenyuan stream salamander, with 40.1% of its transcriptome complete, and the Wyoming toad (Anaxyrus baxteri), at 75% completeness (Carlson et al. 2022), may lack PLCD4 due to incomplete datasets. In contrast, the absence of a PLCD4 transcript in the Cayenne caecilian (Typhlonectes compressicauda), with its 97.6% complete transcriptome, strongly suggests a PLCD4 deletion, but this cannot be definitively verified without a complete genome assembly.

The absence of the PLCζ isozyme in the Xenopus species analyzed in this study poses significant questions regarding PLC-activated calcium release during their fertilization. While a substantial body of evidence has identified PLCZ1 as the primary sperm factor for egg activation (Swann 2022), its absence in Xenopus challenges the universality of this mechanism. Beyond PLCζ, other sperm factors, such as extramitochondrial citrate synthase (CSL or eCS) (Harada et al. 2007, 2011; Kang et al. 2020) and the postacrosomal WWP-domain-binding protein (WBP2NL or PAWP) (Wu et al. 2007), have also been proposed to initiate this process. However, we did not detect a CSL ortholog in Xenopus or any other anuran species. The orthologs we did find contained the targeting sequence characteristic of mitochondrial citrate synthase (CS). By contrast, PAWP is expressed in Xenopus testes (Fig. 3b), but its postulated role in egg activation has been contested in various studies (Nomikos, Sanders et al. 2014; Nomikos, Theodoridou et al. 2014; Nomikos et al. 2015). This leads us to consider alternative candidates like plcd4, which is highly expressed in testes (Fig. 3b). However, sperm from PLCD4 knockout mice induce calcium waves in eggs after intracytoplasmic injection (Fukami et al. 2001).

We finally consider that a soluble sperm protein may not be the sole factor that increases cytoplasmic calcium in the egg. This aligns with the “receptor model” of egg activation, in which a signaling cascade is activated by binding a sperm surface ligand to an egg receptor (Jaffe 1990). Xenopus laevis eggs express a PLCγ isoform (Plcg1), which can be activated by phosphorylation from receptor tyrosine kinases (Sato et al. 2000; Kadamur and Ross 2013), as well as PLCβ isoforms (Plcb1 and Plcb3), which GPCRs activate (Smrcka et al. 1991; Rhee 2001). However, fertilization in X. laevis does not cause observable phosphorylation at the critical tyrosine residue (Y776) of Plcg1, and inhibitors of either PLCγ or PLCβ pathway did not affect the PLC mediated fast block to polyspermy (Komondor et al. 2023). Another alternative explanation is the “calcium bomb” hypothesis, which proposes that a burst of calcium is introduced to the egg with sperm entry, resulting in egg activation and calcium release from the ER (Jaffe 1983). In other model systems, introducing calcium was insufficient to produce calcium release, indicating that other factors are required (Swann and Ozil 1994). Given these findings, it's conceivable that a combination of these mechanisms may be needed to trigger egg activation in Xenopus and other species successfully. Additional studies will be instrumental in identifying the factors required for egg activation, providing new insights into life's earliest events.

Materials and Methods

Ethics Statement

All animal procedures were conducted using acceptable standards of humane animal care and approved by the Institutional Animal Care and Use Committee (IACUC) at the University of Pittsburgh.

Retrieval of Amphibian PLCZ1 Sequences

We searched for amphibian PLCZ1 orthologs in the TSA database using tblastn queries with PLCZ1 protein sequences from the common toad (B. bufo), Plains spadefoot toad (S. bombifrons), and Gaboon caecilian (Geotrypetes seraphini). Additionally, we retrieved the Pleurodeles waltl PLCZ1 ortholog from the TSA hosted by iNewt (Matsunami et al. 2019). To expand our search, we queried amphibian WGS assemblies with tblastn to identify potential PLCZ1 orthologs. For promising hits, we used miniprot (Li 2023) to recover full coding sequences. After pooling sequences encoding potential PLCZ1 orthologs, we validated the candidates through a reciprocal best-hit blastx search against the nonredundant protein database, considering hits with at least 30% query sequence coverage as tentative PLCZ1 orthologs for further analysis.

Phylogeny of Recovered PLCZ1 Orthologs

To validate our initial PLCZ1 classifications, we compared the recovered protein sequences with those of annotated amphibian PLC isozymes, generating a MSA using MAFFT (Katoh and Standley 2013). We processed the alignment output with trimAL (Capella-Gutierrez et al. 2009) to remove poorly aligned regions. For phylogenetic analysis, we used IQ-TREE 2 (Minh et al. 2020), applying the LG + I + G amino acid substitution model (Le et al. 2008) and performing 1000 replicates of UFBoot (Hoang et al. 2018) and SH-aLRT (Guindon et al. 2010) to assess branch support. Unrooted and circular cladogram outputs were produced using Dendroscope (Huson and Scornavacca 2012).

Synteny Analysis of Anuran Genome Assemblies

We retrieved anuran genome assemblies from NCBI for synteny analysis, selecting only those labeled as annotated representative assemblies. Using Genome Data Viewer (Rangwala et al. 2021), we manually inspected these assemblies for conserved synteny blocks around PLCZ1 and PLCD4. Additionally, we obtained unannotated chromosome assemblies to detect PLCZ1 synteny blocks in species from the Pipidae family: the Zaire dwarf clawed frog (H. boettgeri) and the Marsabit clawed frog (Xenopus borealis). Here, synteny blocks were identified using tblastn with protein queries from X. tropicalis (for proteins other than PLCZ1) and S. bombifrons (for PLCZ1). The sequences were annotated with miniprot (Li 2023) and verified through blastx reciprocal best hits. For visualizing local synteny, we aligned sequences in the same gene order, using gene annotations provided by NCBI (for annotated assemblies) or miniprot (for assemblies without annotation).

Animals

Xenopus laevis adults were obtained commercially (Nasco) and kept in a controlled environment with a 12-h light/dark cycle at 20 °C. Sexually mature X. laevis males were euthanized by immersing them in a solution of tricaine-S (3.6 g/L, pH 7.4) for 30 min.

RNA Isolation

To obtain RNA, testes were removed from 3 euthanized X. laevis males and carefully cleaned them to remove fat and vascular tissue. We prepared tissue for RNA isolation by freezing it with liquid nitrogen and grinding it into a powder with a mortar and pestle. We then isolated RNA following instructions provided by RNeasy and QIAshredder kits (QIAgen).

RNA-seq Library Preparation and Data Acquisition

Before preparing the RNA-seq library, we first determined the integrity of the RNA using electrophoresis and measured its concentration by Qubit (Life Technologies). We then used Illumina TruSeq mRNA kit (Illumina) with modified protocol: SuperScript IV (Invitrogen) was used for first strand synthesis, the library was amplified with 10 cycles of PCR, and the amplified library was cleaned up with 35 µL AMPureXP beads (Beckman Coulter). The library was sequenced with 75 bp paired-end mRNA reads on an Illumina NextSeq500 platform with a Mid Output 150 flowcell (Illumina).

RNA-seq Analysis

Sequencing reads were uploaded to the public server at usegalaxy.org (Afgan et al. 2018). The reads were then aligned to the X. laevis genome (version 10.1) using HISAT2 (Galaxy version 2.2.1 + galaxy0) with default settings for paired reads. Next, aligned fragments were mapped and quantified with featureCounts (Galaxy version 2.0.1 + galaxy2) (Liao et al. 2014) using the Xenbase gene model as a reference.

BLAST Search for Unmapped PLCZ1

Command line BLASTx was used to search unmapped RNA-seq reads against a custom database containing PLCZ1 protein sequences. We filtered hits using an E-value cutoff of 10⁻⁶, which produced 360 reads. We then aligned these reads to the translated nucleotide sequence of PLCZ1 from a tree frog P. bahiana, which was excluded from the initial database. This produced a noncontiguous alignment covering residues 170 to 654 of the protein. We then cross-validated the identity of the aligned reads using a multiple query BLAST reciprocal best hits search against amphibian sequences from the nonredundant nucleotide (nt) database.

Comparative Transcriptome Analysis

To compare transcript levels, we took the average of 3 individual experiments and reported transcript abundance as the number of transcripts per million. We deposited this data with the NCBI Gene Expression Omnibus (Edgar et al. 2002), which can be accessed using accession number GSE224304. To ensure that our data could be compared with previously existing data, we matched the transcriptome of our dataset against a GEO dataset for X. laevis (J strain) (GSE73419) (Session et al. 2016). We identified 13,277 genes for direct comparison, which we then ranked and compared using Spearman’s rank correlation in Prism (GraphPad).

To compare the expression of PLC subtypes in the testes of different organisms, we obtained GEO data for X. laevis (J strain) (GSE73419), X. tropicalis (nigerians) (GSM5230669), B. orientalis (GSE163874), and Mus musculus (GSE181426). Datasets were used as provided, except GSE181426, which was annotated and converted to transcripts per million using data from BioMart (19144180). We processed the expression data using a log₂ transformation and created a heatmap using Prism (GraphPad) to visualize the interspecies comparison.

Supplementary Material

Supplementary material is available at Genome Biology and Evolution online.

Acknowledgments

We thank D. Summerville for excellent technical assistance and M.T. Lee for stimulating conversations.

Author Contributions

Rachel E. Bainbridge, Joel C. Rosenbaum, Anne E. Carlson (Conceptualization), Rachel E. Bainbridge, Joel C. Rosenbaum, Paushaly Sau, Anne E. Carlson (Data curation), Rachel E. Bainbridge, Joel C. Rosenbaum (Formal analysis), Rachel E. Bainbridge, Joel C. Rosenbaum, Anne E. Carlson (Investigation), Rachel E. Bainbridge, Joel C. Rosenbaum (Methodology), Rachel E. Bainbridge, Joel C. Rosenbaum, Anne E. Carlson (Project administration), Rachel E. Bainbridge, Joel C. Rosenbaum, Anne E. Carlson (Resources), Rachel E. Bainbridge, Joel C. Rosenbaum, Anne E. Carlson (Visualization), Rachel E. Bainbridge, Joel C. Rosenbaum, Anne E. Carlson (Writing—original draft), and Rachel E. Bainbridge, Joel C. Rosenbaum, Paushaly Sau, Anne E. Carlson (Writing—reviewing & editing).

Funding

This work was supported by National Institute for General Medical Sciences grants 1R35GM153270 and 1R01GM125638 (NIH, USA) to A.E.C. and a grant from the Office of Research, Health Sciences, University of Pittsburgh to J.C.R.

Data Availability

The datasets generated and/or analyzed during the current study are available in the Gene Expression Omnibus, which can be accessed using accession number GSE224304.

Literature Cited

Author notes