All experimental procedures involving model organism A. mexicanum were carried out in accordance with the Czech Law 246/1992 on animal welfare. A. mexicanum animals were from the colony of the Department of Zoology, Faculty of Science, Charles University, Prague, Czech Republic, and all protocols were approved by the Faculty of Sciences of Charles University.

A. mexicanum male and female adults were kept together in an aquarium and after natural stimulation, the females laid eggs. Samples were prepared in two independent experiments and using two different females and males. The gel envelope was first removed from the eggs using tweezers. Eggs were then collected and incubated in sterile 1× Steinberg’s solution containing Pen-Strep (Sigma). Embryos at the 1-, 4-, 64- and 1K-cell stages were embedded in Tissue-Tek O.C.T. Compound (Sakura) with the animal pole oriented at the top. All samples were then stored in the freezer at −80°C.

Samples were subsequently incubated in the cryostat chamber (Leica CM 1950, USA) at −24°C for 10 min and then cut into 30 μm slices along the animal-vegetal axis as shown in Figure 1A. The slices were pooled and equally distributed into 5 tubes. The tube labelling corresponded to embryo orientation: A—extremely animal segment, B—animal, C—central, D—vegetal, E—extremely vegetal.

The samples were homogenized in 300 μl of TRIReagent^® (Sigma-Aldrich, USA) and total RNA was extracted according to the manufacturer’s protocol. LiCl precipitation was performed to remove inhibitors present in the yolky vegetal hemisphere. The concentration of RNA was measured using NanoDrop-2000 (ThermoFisher, USA) and sample quality was assessed using 5200 Fragment Analyzer (Agilent, USA).

The cDNA was prepared using 30 ng of total RNA and RNase-free distilled water (ThermoFisher, USA) in a volume of 5.5 μl and a reaction mixture was added containing 0.5 μl of dNTPs (10 μM each, ThermoFisher, USA), 0.5 μl of oligo-dT and random hexamer (1:1 mixture, 50 μM each, ThermoFisher, USA), and 0.5 μl of RNA-spike (TATAA biocentre, Sweden). The mixture was incubated for 5 min at 65°C and 10 min at 4°C. During the second step, the second mixture was added containing 2 μl of 5xRT Buffer (ThermoFisher, USA), 0.5 μl of RnaseOUT (ThermoFisher, USA) and 0.5 μl of Maxima H Minus Reverse Transcriptase (ThermoFisher, USA). The reaction proceeded as follows—10 min at 25°C, 30 min at 50°C, 5 min at 85°C and cooling to 4°C. The cDNA was diluted to 100 μl using Tris-EDTA buffer solution (Sigma-Aldrich, USA) and stored at −20°C.

qPCR was performed to detect the localization of known genes (list of used primers is attached in Supplementary Table S1). PCR primers were designed using Primer3 (Untergasser et al., 2012). The expected length of qPCR products was 80–120 bp and the annealing temperature was 60°C. Geneious prime (version 2021.2) was used to increase the specificity of designed primers and to avoid targeting RNA isoforms.

qPCR reaction mix with a total volume of 7 μl contained 2 μl of cDNA, 0.29 μl of forward and reverse primers mix (1:1, 10 μM each), 3.5 μl of 2x TATAA SYBR^® GrandMaster^® Mix (TATAA Biocenter, Sweden) and 1.21 μl of RNase-free distilled water (ThermoFisher, USA). qPCR was performed using CFX384 Real-Time System (Bio-Rad, USA) as follows: the initial denaturation for 3 min at 95°C, 45 cycles of denaturation for 15 s at 95°C, annealing at 60°C for 20 s and extension at 72°C for 20 s. qPCR melting curves were analyzed to test the reaction specificity. qPCR data were analyzed using workflow published in Sindelka et al., 2008.

We used 100 ng of a total RNA for library preparation. Ribosomal RNA depletion was performed using Ribocop rRNA Depletion Kit V1.3 (Lexogen, Austria). Libraries were prepared using NEBNext^® Ultra™ II Directional RNA Library Prep Kit for Illumina^® (New England Biolabs, USA). The number of PCR cycles was set at 12 cycles according to the initial RNA concentration. Library concentration was measured with a Qubit 4 Fluorometer (ThermoFisher, USA) and quality was assessed using a 5200 Fragment Analyzer (Agilent, USA). The pooled libraries were sequenced using Illumina NextSeq 500, high-output 150 bp run.

The cDNA of grip2, dnd1, rbpms2 and AMEXTC_0340000004005 was PCR amplified using primers (Supplementary Table S2) designed for a full length cDNA. PCR reaction mix contained 5x Phusion green HF buffer (ThermoFisher, USA), 10 μM MgCl₂ (ThermoFisher, USA), 10 μM dNTP₃ (ThermoFisher, USA), Phusion Hot Start II DNA Polymerase (ThermoFisher, USA), UltraPure dH₂O (Invitrogen, USA) and forward and reverse primers. PCR program run as follows: initial denaturation at 98°C for 30 s, 39 cycles of denaturation for 10 s, annealing at 55°C for 30 s and extension at 72°C for 2 min.

The plasmid pBluescript II KS+ and amplified cDNA were digested using XhoI and NotI (New England Biolabs, USA). 5′-ends of DNA were dephosphorylated using CIP (calf intestinal alkaline phosphatase, New England Biolabs, USA) and then DNA insert was inserted into vector DNA in a ligation reaction. Ligation mix contained T4 DNA Ligase Buffer (New England Biolabs, USA), vector DNA, insert DNA, RNase-free distilled water (ThermoFisher) and T4 DNA ligase (New England Biolabs, USA). The mix was incubated overnight at 16°C. Then, the reaction was stopped at 65°C for 10 min.

NEB 5-alpha competent E. coli was used for a transformation (High efficiency, # C2987I, New England Biolabs, USA) according to the manufacturer’s protocol: High Efficiency Transformation Protocol (C2987H/C2987I). After that, the individual clones were sequenced to screen for the presence of the expected sequences. Correct plasmids were purified using Plasmid Midi Kit (Qiagen, Germany).

Plasmids were linearized in the restriction digest reaction. The mix contained 7 μg of plasmid DNA, NEB restriction enzyme (New England Biolabs, USA) and 10x NEB buffer 3.1 (New England Biolabs, USA). Mix was incubated overnight in a 37°C water bath. Linear DNA was cleaned up using QIAquick PCR Purification Kit (Qiagen, Germany). During the transcription reaction, we mixed 2 μg of linear template, 4 μl of 5x transcription buffer (Agilent, USA), 2 μl DIG RNA labelling mix (Roche, Switzerland), 2 μl of the T7 polymerase (Agilent, USA) and water in total volume of 20 μl. The mixture was incubated in a 37°C water bath for 3 h. Then, the mix was cleaned using LiCl. Finally, the RNA probe quality was tested using formaldehyde gel.

In situ hybridization was performed on whole mounts as described previously (Soukup et al., 2021) with slight modifications. Briefly, rehydrated A. mexicanum albino (d/d) embryos were digested in 60 μg/ml Proteinase K in PBS, fixed in 4% formaldehyde + 0.2% glutaraldehyde for 30–120 min, transferred into hybridization solution (50% formamide, 1x Denhardt’s, 1 mg/ml yeast RNA, 0.1% Tween-20, 10% dextran sulfate, 1x salt solution containing 0.2 M NaCl, 8.9 mM Tris-HCl, 1.1 mM Tris base, 5 mM NaH₂PO₄.H₂O, 5 mM Na₂HPO₄ and 5 mM EDTA), and incubated overnight in hybridization solution containing RNA probe (1:1,000–1:100). Next day, the specimens were washed several times in post-hybridization solution (50% formamide, 4 × SSC, 0.5% Tween-20) and transferred via MABT buffer (100 mM maleic acid, 150 mM NaCl, 0.1% Tween-20) into blocking solution (2% blocking reagent, 20% sheep serum, in MABT buffer). Following blocking, the specimens were incubated overnight in the blocking solution containing alkaline phosphatase-conjugated antibody against DIG (Roche, 1:3,000) at 4°C. The specimens were washed several times in the MABT buffer. Following the overnight MABT wash, the samples were transferred into NTMT buffer (0.1 M Tris, 0.1 M NaCl, 0.05 M MgCl₂, 0.1% Tween-20) and incubated in BM Purple substrate (Roche) at 4°C until desired signal developed.

RNA-seq reads were processed as previously described in Naraine et al. (2022). Adaptor sequences and low quality reads were removed using TrimmomaticPE (v. 0.36) (Bolger et al., 2014) using the parameters, “HEADCROP:12 ILLUMINACLIP:∼/TruSeq-PE3.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36.” Mitochondrial RNA reads (GenBank id: AY659991.1) and any remaining rRNA reads were removed using SortMeRNA (v. 2.1b) (Kopylova et al., 2012). The reads were then pseudo-aligned to the A. mexicanum transcriptome AmexT_v34 (Nowoshilow and Tanaka, 2020) using kallisto (v. 0.43.1) (Bray et al., 2016). The data were deposited in the National Center for Biotechnology Information’s Gene Expression Omnibus (GEO: GSE240796).

Raw counts were initially filtered to keep transcripts with counts greater than 30 in at least one sample. DESeq2 (v. 1.32.0) (Love et al., 2014) was used to normalize the counts using the median-of-ratios method followed by differential expression analysis to determine differential localization of transcripts along the animal-vegetal sections in the 1-cell, 4-cell, 64-cell, 1K-cell stages; changes in the sectional profile across all stages; and changes in the total transcript across the stages. The median-of-ratios normalization method was used to focus primarily on the extremely localized transcripts.

The following DESeq2 design models were used:

1) Alteration between the sections at the same stage:

The Principal Component Analysis (PCA) of the top 500 variable transcripts was assessed for the presence of any outlier samples. Differentially localized transcripts (DLTs) were defined as those with an adjusted p-value (padj) value less than 0.01 and also a total transcript count greater than 20 within at least one stage. Human gene symbols were assigned based on the previous ortholog analysis from Naraine et al., 2022, whereby A. mexicanum gene symbols were either matched against all known Homo sapiens gene symbols or derived from the similarity between its protein sequences as compared against the H. sapiens proteome using the reciprocal best alignment heuristic tool Proteinortho (v. 6.0.9) (Lechner et al., 2011).

The spatial expression profiles of the DLTs were then characterized into five discrete categories (extreme animal, animal, central, vegetal, extreme vegetal) based on the parameters previously described in Naraine et al., 2022. The central category represents a new category of localized maternal transcripts that was observed primarily in the eggs of A. mexicanum and A. ruthenus (Naraine et al., 2022). The other four profile parameters remained unchanged from the previous publication describing sub-location in X. laevis (Sindelka et al., 2018). DLTs that did not fit into these defined parameters were labelled as unclassified. Profile changes in the DLTs across the whole stages and within the stages were assessed using the degPatterns function from DEGreport (v. 1.28.0) package (Pantano et al., 2023). Profiles where the DLTs showed a fold change of 3x or 2x difference between either the stages or the sections respectively were selected. The validity of the profiles was verified using optCluster (v. 1.3.0) with the “Diana” clustering algorithm (Sekula et al., 2017).

Gene ontology terms associated with the genes were obtained using online software g:Profiler (access date: 14/02/23) using the default parameters of the annotated human reference, multiple testing correction using g:SCS threshold with a cutoff of 0.05 (Raudvere et al., 2019). Gene ontology terms clustering and removal of redundant terms were done using Revigo (access date: 14/02/23) with the default parameters of the whole UniProt reference database and SimRel semantic similarity measurement (Supek et al., 2011).

We analyzed the presence of motifs unique within the UTRs of A. mexicanum PGC markers that might be responsible for its degradation. This was done by comparing their 3′ and 5′UTRs against shuffled versions of their sequences or the UTRs from X. laevis PGC genes. We also analyzed for enrichment of motifs within the vegetal and animal transcripts using the previously observed vegetal and animal motifs that were detected in A. mexicanum oocytes (previously published by Naraine et al., 2022). Analyzing the 3′UTR sequences, we also checked for de novo motifs that might be enriched selectively within the degraded transcripts versus the de novo transcripts and also between the animal and vegetal groups.

Motifs were detected using the STREME software (v. 5.5.2) (Bailey, 2021) under the following conditions: p-value < 0.05 and motif width = 6 to 25. The motif enrichment in primary sequences compared with control sequences was assessed using AME (e-value ≤ 0.05) (v. 5.5.2) (McLeay and Bailey, 2010). FIMO software (p-value ≤ 0.0001) (v. 5.5.2) (Grant et al., 2011) was used to scan the identified motifs against primary and control sequences. Using FIMO we obtained information about motif position within the sequences, and it was used to assemble a map of the motif distribution.

The enrichment of the previously observed vegetal and animal motifs from the paper (Naraine et al., 2022) was analysed against the UTRs of the animal and vegetal transcripts observed during embryogenesis. Significantly enriched motifs were deemed as those that gave an AME e-value ≤ 0.05 and a 3x fold enrichment in the vegetal UTRs versus the animal UTRs. In the case of the occurrence of the same or very similar motifs, we selected the motif with the lowest e-value.

Using only the PGC dataset, we continued with the analysis of putative regulatory elements in the UTRs that might explain its temporal degradation. To identify known protein binding sites and RNA binding proteins (RBPs) within the 3′UTR, we used Scan For Motif (access date: 16/03/23) with the datasets from TransTerm (E-value <= 0.175 per thousand bases) and RBPDB (E-value <= 0.001) and selected results with E-value thresholds < 0.001 (Biswas and Brown, 2014). The motifs recognized by these RBPs were compared against the list of de novo motifs detected by STREME using the comparison tool Tomtom (v. 5.5.2) (Gupta et al., 2007). The position probability matrices for the motif of these RBPs were downloaded from several online databases: CIS-BP-RNA (v. 0.6) (Ray et al., 2013), RBPDB (v. 1.3.1) (Cook et al., 2011) and oRNAment (access date: 12/06/23) (Benoit Bouvrette et al., 2020). Within 3′ and 5′UTR, the BEAM software (v. 1.6.1.) (Pietrosanto et al., 2018) was used to find RNA secondary structures using a p-value threshold cutoff of 0.01. The BEAGLE software (access date: 23/05/23) (Mattei et al., 2015) was then used to identify a conserved secondary structure either shared across all the A. mexicanum PGCs or unique to A. mexicanum relative to the PGCs in X. laevis, using a p-value threshold cutoff of 0.01. BRIO (p-value < 0.05) (Guarracino et al., 2021) was used to identify known sequences and structure RNA-binding motifs that are recognized by RBPs in the UTRs of H. sapiens and M. musculus from PAR-CLIP, eCLIP and HITS experiments. The AURA database (v. 2.7) was then used to screen the identified RBPs for their selective preference to either the 3′ or 5′ UTR regions (Dassi et al., 2014). For the discovery of putative binding sites for miRNA within the 3′UTR of the PGC markers, miRDB (v. 6.0) (Chen and Wang, 2020) was used with the human dataset as a reference.

To describe RNA localization during early embryogenesis of A. mexicanum, we collected the stage of a fertilized egg (1-cell stage), 4-, 64- and 1K-cell (early blastula), which are expected to be pre-embryonic genome activation (MBT).

Embryos were sectioned along the animal-vegetal axis using the TOMO-seq method (Figure 1A). The results were first analyzed using PCA of the 500 most variable transcripts (Figure 1B). PC1 showed high variability during early development, mainly between the 1-cell stage and late embryos. PC2 revealed the clear distinction among individual sections in all stages. The variability in the sections decreased as the embryos progressed towards the later stages, indicating a disruption or reduction in the original asymmetrical gradients. Additionally the largest sectional differences were observed between the D and E sections, while minimal between the A and B sections.

The average number of coding genes identified at each stage was 28909 (using threshold >30 transcript counts in any sample). The complete dataset containing genes with asymmetrical distribution are listed in Supplementary Table S3. More than 2,200 DLTs (padj < 0.01, using threshold >20 transcripts in at least one stage) with sectional changes per stage were identified in embryos at 4, 64 and 1K-cell stages and twice as many were identified at the 1-cell stage (4,076 DLTs). The diagram in Figure 1C shows the overlap of the shared DLTs across the stages with 1,216 shared in all analyzed stages. Most of the DLTs were classified into one of the five localization categories: extremely animal, animal, central, vegetal and extremely vegetal (Fig. 1Da). In extremely animal and animal categories, we identified 284–675 DLTs in each stage (Fig. 1Db). Majority of DLTs are localized in vegetal or extremely vegetal sections at the 1-cell stage—1,114 and 1,083. The number of extreme vegetal DLTs dramatically decreased during development and only 240 were found at the 1K-cell stage.

We performed RT-qPCR validation of a few members within the extremely vegetal and animal transcript categories. In the extremely vegetal category, we confirmed the localization profile of grip2 and dnd1 (Supplementary Figure S1A), and in the animal category, we confirmed the localization of ankhd1 and akt2 (Supplementary Figure S1B).

In total 4,850 DLTs were observed to change their sectional profile or its sectional amplitude across the analyzed stages. 781 of these DLTs had a ≥2x fold change between a given section and showed two profile alteration processes—alteration of vegetal profiles and alteration of animal profiles (Supplementary Table S4).

Alteration in transcript count caused by synthesis or degradation was detected for 6,811 DLTs (padj <0.01, >20 transcripts in at least one stage). Out of these, we observed 960 DLTs with at least 3-fold change across individual stages (Supplementary Table S6). Two main groups were created reflecting de novo synthesis and degradation.

Transcripts degradation detected using the TOMO-seq approach was validated by RT-qPCR for 2 DLTs—plin2 and velo1 (Supplementary Figure S1C). De novo synthesis of mok and nynrin-like was also validated by RT-qPCR (Supplementary Figure S1D).

For a deeper understanding of the PGC markers degradation process, we performed motif analysis searching for any conserved regulatory sequence within the 3′ and 5’ UTR sequences. Using data from previous publication (Naraine et al., 2022) and Xenbase.org we obtained UTRs of 5 transcripts in each organism (Listed in Supplementary Figure S2A). There were no statistically significant motifs enriched within the PGCs RNAs of the A. mexicanum versus those from the X. laevis. However, this is most likely due to the low number of genes used for the analysis (∼5 genes). Due to this limitation, FIMO was instead used to assess whether the detected motifs were found exclusively or in a high proportion within the UTRs of the A. mexicanum PGCs versus the other models.

The de novo motif analysis within the 3′UTR of PGC marker genes revealed 7 motifs enriched in the A. mexicanum compared with either the X. laevis marker genes or the shuffled sequences (Figure 6A). Most of the motifs were uridine rich. 5 of the identified motifs were present exclusively in A. mexicanum but not in X. laevis. Then, the motif presence was also estimated in the other model organisms—M. musculus, H. sapiens, D. melanogaster, A. ruthenus and D. rerio—and only 1 motif (Motif 6) was exclusive to A. mexicanum. We observed no significant enrichment of these 7 motifs when comparing the de novo transcripts against the degraded ones, or the vegetal sectional changes against the animal sectional changes.

In total, we identified 10 motifs conserved within the 5′UTR of PGC markers in A. mexicanum and none of them were present in X. laevis or the shuffled sequences (Figure 6B). Most of the motifs were enriched with either cytosine or guanosine. Two motifs (Motif 1 and 6)—had CAC core. All the motifs were scanned against PGC sequences of the 4 model organisms and 4 motifs were unique exclusively to A. mexicanum. We observed no significant enrichment of these 10 motifs when comparing the de novo transcripts against the degraded ones, or the vegetal sectional changes against the animal sectional changes.

We used the Scan for motif software to identify putative RBP sites within the 3′UTR of PGC markers. The analysis revealed several RBPs which can affect RNA stability, degradation and translation, such as PUM2, KHSRP and ZFP36. However, all the RBPs affecting RNA stability were also detected in X. laevis PGC markers (Complete list of RBPs in Supplementary Table S8). To compare de novo motifs with motifs recognized by known RBPs we used the Tomtom tool, which helped us to identify 2 motifs which resembled previously discovered motifs. Motif 5 was similar to RBFOX1 (UGCAUG) binding sites and motif 6 was similar to PUM2 (UGUA) and ZFP36 (UUAUUUAWK) (Figure 6A).

BRIO software was used to search for known primary sequences and secondary structures associated with RBP binding motifs within the 3′ and 5′ UTRs of A. mexicanum PGC markers. We identified 5 (3′UTR) and 17 (5′UTR) RBPs binding motifs unique to A. mexicanum compared with X. laevis (Supplementary Table S9). Using the available resources, we selected 2 (3′UTR) and 9 (5′UTR) proteins potentially playing a role in transcript degradation (Figure 6C). BEAM was used to find structural motifs present in the RNA secondary structure, but no motif was shared in all 5 A. mexicanum PGC markers within 3′ or 5′UTR.

As it is known, miRNA can induce mRNA degradation or translational repression (reviewed in O’Brien et al., 2018). Therefore, we searched for miRNAs that may target the A. mexicanum PGC transcripts using the online miRNA database miRDB. We discovered 41 miRNAs present only in A. mexicanum, but none of them was shared in all PGC genes. Similarly, 50 miRNAs were specific only for X. laevis but were not shared in all its PGC genes (Supplementary Table S10).

Analysis was carried out to determine if previously conserved animal and vegetal motifs detected in the oocyte of the A. mexicanum (Naraine et al., 2022) can also be seen preferentially enriched within the vegetal and animal transcripts during embryogenesis. Motifs were scanned using FIMO against the following datasets of developing A. mexicanum embryos: animal, vegetal, sectional profile altering, stage altering de novo and degraded DLTs.

Within the 3′ and 5’ UTR of transcripts from the A. mexicanum embryos, we observed the enrichment of the previously detected putative localization motifs that were found shared between both A. mexicanum and X. laevis oocytes and also those that were unique only for A. mexicanum oocytes (Supplementary Figure S3) (Naraine et al., 2022). The conserved motifs were detected only in the vegetally localized DLTs of A. mexicanum embryos and were mostly enriched with CAC or guanine-rich sequences.

The map of identified de novo motifs, miRNAs, RBPs and motifs from previous publication (Naraine et al., 2022) was assembled for selected PGC markers (Supplementary Figure S4).