Basal Xenobot transcriptomics reveals changes and novel control modality in cells freed from organismal influence

In this study, we ask the question, how does the transcriptome change in response to cellular collectives acquiring non-evolved embodiments and behaviors in absence of any genomic editing, drug treatments, or heterologous interventions—to understand the top-down influence of morphology over transcriptome. Our goal was to utilize a synthetic living system to identify unique transcriptomic changes in it, in comparison to its native (wild-type) context. For this, we exploited the Xenobots platform, which includes multiple different kinds of autonomously-moving biobots all derived from Xenopus embryonic cells32–35,37,107. We analyzed the basal Xenobots, which are autonomously-motile, self-assembling multicellular entities derived from Xenopus ectodermal explants (animal caps)33. The Xenopus ectodermal explants have a rich history of serving as a tractable model system for investigating epidermal cell fate specification and cell fate plasticity in response to various chemical stimuli49–55. In addition, mucociliary organoids derived from them have helped understand mucociliary epidermis dynamics56–61.

Recent transcriptomic analysis has looked at epidermal cell fate specification during transformation of these ectodermal explants to the autonomously-moving entities (basal Xenobots)56,62. However, in this study, we investigate the transcriptome of these Xenobots in comparison to their native wild-type embryonic tissue using a different transcriptomic analysis beyond epidermal cell fate specification, to understand how these Xenobots’ morphology and behavior affect their transcriptomic changes. We found significantly greater inter-individual gene variability in these Xenobots, suggesting adaptive exploration of the transcriptional space to adapt to their novel embodiment and associated stressors. Gene ontology and network analysis of unique (non-epidermal marker) transcriptomic changes showed both expected shifts in multicellular organization and immune/stress response along with surprising changes in metabolism and sensory perception of sound and mechanical stimuli. Indeed, these autonomously motile basal Xenobots respond to acoustic stimulation by changing their motion behavior while motile age-matched embryos and passive non-motile Xenobots do not. Physlostratigraphic analysis showed transcriptional changes reflecting evolutionary ancient genes. We also found significant enrichment of thanatotranscriptomic genes, suggesting signatures of organismal death in them. Overall, all these transcriptomic changes in these Xenobots suggest that alteration of morphology and nascent emergent life history were enough to strongly affect the transcriptome (morphology/behavior driving gene expression).

Since these Xenobots are raised in different media concentration (0.75X MMR) as opposed to embryos (0.1X MMR), we first wanted to see if there are any major transcriptome changes in response to the different media concentration. Since the ectodermal explants were excised from embryos at stage 9 and incubated in 0.75X MMR for making these Xenobots, we incubated embryos in 0.75X MMR from stage 9 onwards until age-matched stage to mature day 7 Xenobots. Raising embryos in 0.75X MMR after stage 9 remarkably showed normal embryonic development with no developmental morphology defects, and embryos looking similar to their sibling raised in 0.1X MMR. Comparing transcriptome of embryos raised in 0.1X MMR to 0.75X MMR showed no major transcriptome differences (only 6 transcripts out of ~24,500 transcripts showed significant differences) (Supplementary Fig. 1). This was very surprising and suggests that the embryonic transcriptional regulation is very well shielded from external environmental changes. This might be due to the presence of tight junctions in the outmost ectodermal layer, which compartmentalizes the embryos from their external environment. If embryos are raised in 0.75X MMR from pre-stage 9, then they do not undergo normal gastrulation and do show developmental defects, suggesting as previously documented, that the pre-gastrulation embryos are much more susceptible to environmental differences than a post-gastrulation embryo.

To deal with false positives, we used scRNAseq database of Klein tools79. However, the caveat is that we had to use scRNA data from stage 22 embryos, whereas our age-matched embryos were stage 35. This stage 22 scRNA dataset is the best accessible results we have as there are no stage 35 scRNAseq datasets which we could use to remove epidermal progenitor transcripts. However, by stage 22, epidermal cell differentiation has already occurred62,79. In addition, to eliminate false positives, we took an extremely stringent approach by removing all genes expressed in epidermal progenitor cells, multiciliated cells, alpha and beta ionocytes, and goblet cells from our Xenobot dataset. In fact, this stringent curation may have resulted in the loss of some true positives, but it was worth it to give us good confidence in our results. Overall, our highly stringent list of Xenobot transcripts is a good starting point to look into unique transcriptional aspects of these Xenobots and to begin understanding how genome/transcriptome relates to novel morphologies. In future, single-cell RNA sequencing of Xenobots might allow us to look at regional and cell-type-specific expression patterns in Xenobots and how ectodermal progenitor cells distribution and profile diverge between Xenobots and embryos.

We compared and analyzed the transcriptome of these Xenobots with their age-matched Xenopus embryos to determine effects of altered morphology and nascent emergent life history on transcriptome. Our results give an indication of the kinds of transcriptional changes induced by change of morphology and behavior (no genomic editing or transgenes). We found at least 537 transcripts uniquely upregulated in these Xenobots (Supplementary Data 4). These transcripts indicated enriched functions such as locomotion apparatus assemblies (cilia-based motion and motor proteins), structure formation and multicellularity, stress and immune response, shift in metabolism, and sensory perception of sound and mechanical stimuli (Fig. 3, Supplementary Fig. 3, and Supplementary Data 5 and 6). The majority of these shifts were towards evolutionary older systems, which likely get suppressed when they are part of the organism. Phylostratigraphic analysis shows the transcripts uniquely upregulated in these Xenobots compared to the wild-type tissue are enriched for evolutionarily ancient genes and systems (Fig. 5). This is fascinating and suggests that cellular functional plasticity can roll backwards along the “ontogeny recapitulates phylogeny” axis, utilizing gene ensembles and systems from their distant evolutionary past that might be relevant in their new configuration and environment. Future work will examine possible links to similar phenomena (a shift toward ancient, unicellular transcripts) in cancer108 and potential implications for the relationship between cancer and multicellularity109,110.

The transcripts uniquely upregulated in these Xenobots compared to the wild-type tissue also showed enrichment for thanatotrascriptomic genes71,73,94–97,111, which are transcripts awakened in cells and organs after the overall organism dies. These transcripts mainly included immune/stress response and insulin signals. What this suggests is that these Xenobots have some sort of memory/information of them not being part of “living organism”. One possibility is that this might be evidence that the cells experienced the Xenobot construction event as the death of their parent organism; it is also likely that the thanatotranscriptome is limited by the fact they retained multicellularity; perhaps dissociating the cells and having them reboot multicellularity (as occurs in a different type of Xenobot33) might show a much stronger thanatotrascriptomic response. One challenge in this field has been that it was not possible to test whether these death-induced genes were functionally active. Indeed, what kind of function could even be assayed at the end of life? The Xenobots platform, as a whole, offers a path forward, because in future work, the thanatotranscriptome could be down-regulated via RNAi or morpholinos, to examine whether this would induce any negative effects on their continued survival or Xenobot functionality.

When cells are challenged and need to make quick changes to their gene expression to adapt to novel conditions, they exhibit exploratory adaptation by changing many different genes simultaneously, and different cells change different sets of genes to overcome the same challenge66–69. When comparing the overall transcriptome of the basal Xenobots and age-matched Xenopus embryos using PCA, the three Xenobots samples showed much wider separation compared to the three age-matched Xenopus embryo samples (Supplementary Fig. 2). Further, our analysis of gene count variability showed significantly higher inter-individual gene variability in these Xenobots compared to age-matched Xenopus embryos (Fig. 2). This increase in inter-individual gene variability in these Xenobots suggests that these Xenobots might be exploring the transcriptomic space (of possible gene expression profiles) to discover transcriptomes that allow them to improvise and adapt to their new embodiment. This is supported by the top ten most variable transcripts in these Xenobots which show prominent presence (three genes out of ten) of oncomodulin that encodes a protein usually found only during early embryonic cells and tumor cells and performs widely distinct functions78 including maintaining sensory perception, tissue regeneration, immune response, and strong antioxidant properties78. Other genes in the top ten most variable transcripts in these Xenobots include: tenomodulin which is involved in tissue maintenance by supporting stem cell renewal and preventing senescence and aging112–114, adhesion proteins neurotrimin115 and salivary glue protein-3 like protein116, A.

superbus venom factor 1 an immune regulator117, somatostatin receptor 4 a regulator of hormone and secretory protein secretion118, chloride intracellular channel 4 that belongs to an extremely conserved (across prokaryotes and eukaryotes) family that can perform different and independent functions at the membrane and the cytoplasm including fundamental cellular processes such as regulating mitochondrial function, exosome communication, membrane trafficking, and pH maintenance119,120, and tyrosine aminotransferase which is involved in gluconeogenesis— an alternate metabolic pathway121.

All of these functions are likely to be useful in the capacity to adjust to novel embodiments and challenges. Xenobots are an intriguing model in which to characterize problem-solving during traversal of transcriptional, physiological, and other unconventional spaces88,122. Individuals may explore the option space for many possible paths to adaptations to new ways of existence. Estimation of total cell number of a cohort of Xenobots reveals that they are remarkably similar in total cell number (Supplementary Fig. 4), ruling out large variability in cell numbers as a reason for transcriptomic differences (as also do prior studies showing embryos’ remarkable invariance to cell number123,124). Alternatively, the greater transcriptomic differences among these Xenobots may be a reflection of heterogeneity introduced to these Xenobot populations by transcriptional memory of the unique “lived experiences” of each Xenobot125–127. Finally, it has been observed that embryos in groups assist the morphogenesis of others and develop more uniformly with fewer defects128. While the embryo groups in the present study were small (10 individuals), they may still have some capacity to synchronize their transcriptomes that these Xenobots lack. In future studies, a spatial transcriptomic analysis of individual Xenobots would provide a better understating of this adaptive exploration process.

We saw not a single transcript overlap between the uniquely upregulated transcripts in these Xenobots in comparison to wild-type tissue and genes representing the mesoderm, endoderm, and axis patterning clusters (Supplementary Data 7), suggesting that there was no “contamination” of the explanted ectodermal tissue from which these Xenobots were derived with mesoderm, endoderm or axis patterning genes. Likewise, it did not support transdifferentiation of ectoderm cells into other lineages.

Overall, both functional enrichment analysis and the network clustering analysis showed largely overlapping functional clusters (Fig. 3, Supplementary Fig. 3, Supplementary Data 5, and 6). Even with the stringent subtraction conditions, including subtraction of all multiciliate cell genes, we still saw cilia assembly and motor proteins as major enrichment categories in both gene ontology analysis and network clustering analysis. This suggests that there might be significant upregulation of generation/making of multiciliate cells in these Xenobots compared to embryos. This is supported by our previous observation that Xenobots have >2 fold higher density of multicilated cells compared to age-matched embryos33. We did not observe much difference in cilia length and cilia number between Xenobots and age-matched embryos (Supplementary Fig. 5). In addition, the multiciliate cells in embryos may be much more stable and long-lasting compared to those Xenobots, which might be making them at a higher rate due to rapid turnover. This is partly supported by the observation of both proliferation and cell death categories being upregulated in these Xenobots, suggesting a higher turnover of cells. It is also possible that there are some unique and subtle differences in cilia assembly and adjoining motor proteins and their functions in these Xenobots’ morphology configuration. Overall, there is definitely a major shift in cilia assembly and motor proteins and perhaps multiciliate cell dynamics in these Xenobots; this remains an area of active investigation.

Subcategories involved in tissue building, repair, and multicellular organization were ATP binding, calcium binding, and subcategories in stress and immune response signals, mainly included cytokines, heat-shock proteins, and MMPs. It has been postulated that stress/immune response is important indicator of deviation from homeostasis and target morphology states and involved in moving the system by repair and tissue building processes to reach its final state88–92. It was interesting to see V-ATPase among these genes, as it is not only deeply conserved (evolutionarily ancient) but has also been postulated to be one of the key elements required for generation of multicellularity129–131, morphogenesis132–135, and cancer136–138. These genes represent candidates for control knobs that will be tested for their ability to enable guidance of the self-organization process to produce other types of Xenobots of different shapes and morphologies.

It was surprising to see sensory perception of sound and mechanical stimuli as an upregulated category in these Xenobots. Initially, we thought it might be erroneous or cross-assignment of cilia-related transcripts, but transcripts in this category were not cilia-related. They indeed contained transcripts known to be involved in sound and mechanical stimuli perception such as: GJB2 which is a gap junction protein involved in sensory perception (among other functions) and mutational loss of this protein leads to hearing loss in humans and mice139–141, stereocilin142—the mechanoreceptor for sound waves/vibrations, proteins involved in inner ear membrane formation, inner ear ion channels, and many other stereocilia (inner ear hair cells) specific proteins. Importantly, the expression data were experimentally confirmed in a functional experiment. We tested the ability of these basal Xenobots to respond to acoustic vibration stimulus (Fig. 4 and Supplementary Movies 1–6). Observable changes in their motion behavior (Fig. 4 and Supplementary Movies 5 and 6) confirm that they do indeed respond to acoustic stimulus. This change in motion behavior is not due to physical movement of medium or simple mechanical effect of vibration on cilia function (Fig. 4 and Supplementary Movies 1–4), as indicated by controls (immotile bots and ciliated age-matched normal Xenopus embryos, neither of which responded to the stimulus).

Future work will detail the specific mechanisms leading from stimulus to change in behavior and characterize additional motility patterns that may be inducible by diverse stimuli. Importantly, the field of morphological computing reveals how the physical properties of bodies serve as an integral part of their cognition143. The Xenobots’ rich mix of ciliary action (itself a highly active multiscale system), internal bioelectric and calcium physiology, and fluid dynamics makes Xenobots an exciting model system in which to unravel morphological computing and behavioral responses in novel morphologies with wild-type genomes144–149.

Xenobots' behavioral response to acoustic stimulus suggests several points for consideration. First, these basal Xenobots express high levels of acoustic-related genes and are sensitive to such stimuli, whereas age-matched Xenopus embryos are not. The induction of whole cluster of transcriptomic machinery sufficient for a functional behavioral response, in a wild-type genome (in the absence of synthetic circuits or added induction by novel cell populations), may shed light on evolutionary developmental biology consequences of changes of morphology and nascent emergent life history. This in turn suggests a future research program to attempt to identify novel sensory-behavioral capabilities in other minimal, organoid, and bioengineered living constructs, which may heretofore have gone unnoticed. Future studies in organoid behavior and their potentially novel transcriptomics can offer new ways to exchange information with them, as shown here with Xenobots using acoustic stimuli identified through transcriptomics. Similar implications also apply to patient-derived-xenografts150. Second, from the perspective of biorobotics and the use of this model system as useful synthetic living machines, the ability to reversibly modulate their behavior via acoustic stimuli opens fascinating possibilities for control. Work is ongoing to map the frequency, waveform, and amplitude domains and link them to predictable changes in behavior (both via motility and via gene expression and physiological state). Real-time modulation of the hydroacoustic parameters in closed-loop controllers is a possible path forward for guided functionality in bio-robotics in Xenobots and beyond.

What would the role of such sensory perception of sound and mechanical stimuli be in such an aneural system like these Xenobots? And how does that regulate its morphology and movement behavior? A role of acoustics in aneural systems is not unprecedented. Recent studies have shown that coral larvae, which are aneural masses of cells with cilia-based movement similar to these Xenobots, perceive reef sounds (distinguishing healthy vs unhealthy reefs and mother reefs vs other healthy reefs), in order to migrate preferentially151–154. It is not known how (mechanism of action) they are able to achieve this acoustic perception, and it remains an area of investigation. Thus, it is possible that such acoustics-based sensing and movement is one of the evolutionarily ancient mechanisms of perceiving the environment and reacting to it via motion before nervous system evolved, leading to coordinated movement behavior. Given the presence of machinery for perception of sound and mechanical stimuli in these Xenobots, acoustic perception and response similar to coral larvae should be examined in Xenobots.

Although it is well known that yolk platelets provide the energy source and essential nutrients for Xenopus during embryonic development155–157, little is known about the metabolism happening within the embryos. Some evidence suggests that during early Xenopus embryonic development (pre-gastrulation), the pentose cycle and glutamate-aspartate cycle are predominant158–160. The Embden-Meyerhof pathway and Krebs cycle, although functional, contribute little during pre-gastrulation stages, but their relative contribution shows major increase after gastrulation158–160. Interestingly, in the basal Xenobots there seems to be upregulation of different metabolic processes, such as gluconeogenesis (which has been shown not to occur in embryos159) and ketone metabolism and metabolic process breaking down cholesterol and steroids. These were associated with presence of many liver-related transcripts involved in metabolism, although these Xenobots have no liver. It is interesting to see upregulation of gluconeogenesis and ketone metabolism in these Xenobots which are highly conserved metabolic processes of generating glucose from non-carbohydrate carbon substrates and is found in microorganisms, bacteria, fungi, plants, and animals161–167, suggesting that metabolically these Xenobots might be shifting to an ancient mechanism to meet their energy requirements. This is perhaps due to the majority of yolk (energy source of embryos) being on the vegetal side of the embryo, and since these Xenobots are derived from the animal side, they may have disproportionately lower yolk, thus needing to manage energy use given the limited amount of energy source.

Interestingly, V-ATPase, the deeply conserved proton pump known to be involved in many functions (generating multicellularity, morphogenesis and cancer129–138) is also known to regulate mitochondrial metabolism168,169 and was found to be upregulated in these Xenobots as part of metabolism cluster, thus suggesting an overall shift in the metabolism in these Xenobots. Further detailed and deep investigation first into the metabolism and energy usage of embryos and then of these Xenobots will help reveal the shifts in metabolism. Understanding this will allow us to manage the energy requirements and control the overall nascent emergent life history of these Xenobots. Along these lines, it has been shown that external supply of glucose can drastically increase the life-span of these synthetic biobots33.

From the gene ontology and network analysis categories, it seemed that the basal Xenobots had activated evolutionarily ancient functions. To test the hypothesis that these Xenobots were acquiring more ancient transcriptomes as a result of their morphogenesis and nascent emergent life history, we performed a phylostratigraphic analysis. Here we compared the control transcripts (all transcripts from epidermal progenitor cells + alpha and beta ionocytes+, multiciliated cells + goblet cells) with transcripts uniquely upregulated in these Xenobots. The controls and Xenobot transcripts were equally enriched until Vertebrata strata. However, although the controls have a much larger transcript dataset (3374), the Xenobot transcripts (1450 uniquely upregulated transcripts) were much more enriched than controls in the older strata Bilateria, Eumetazoa, and Metazoa. This again suggests that these Xenobots have upregulated transcripts and programs that are evolutionarily ancient. Along those lines, the cellular physiology of these Xenobots appears to be closer to evolutionary ancient organisms with a single layer of differentiated cells with no nervous system and vasculature, and cilia-based movements. The phylostratigraphic analysis also shows that these Xenobot genes are evolutionarily basal and found in organisms with simpler body plan. This fits in line with the observation of shifts in metabolism, sensory perception process, and multicellular organization towards evolutionary ancient processes. This indicates that when cells and tissues are liberated from their organisms, they can tap into their evolutionary history to kickstart ancient processes and adjust their energy requirements to reach novel morphology goals.

Taken together, the results of this study reveal a control arc that is complementary to the conventional control of morphology by gene expression: to some extent, the genome is being used as a “resource book” by the cellular collectives to kick start various processes as per the conditions and needs of atypical multicellular configurations. Plasticity in form, function, and molecular physiology is gaining interest, as it impacts on basic evolutionary biology, biomedical applications, and synthetic bioengineering89,170–176, in way that augments current bottom-up approaches of DNA engineering and synthetic biology circuits. We see this work as part of a future roadmap in which a wide variety of synthetic systems are examined to see ways in which living organisms adapt to new environments on-the-fly and solve evolutionarily novel problems in transcriptional, morphological, and physiological spaces. Discovering how groups of cells manage their molecular-biological and metabolic resources to achieve adaptive ends in novel scenarios will facilitate the ability to control and generate different morphologies from the same genome by reshaping the morphological landscape to create new attractor states, control their behavior and responses, and manage their nascent life histories. Such knowledge will serve biomedical purposes (controlling multicellularity during regeneration, birth defects, and cancer), advance biorobotics and bio-AI (design systems with morphology and behavior to solve specific problems), and explore the space of possibilities of life-as-it-can-be.

All experiments were approved by the University Institutional Animal Care and Use Committee (IACUC) under the protocol number M2023-18. We have complied with all relevant ethical regulations for animal use. Xenopus Laevis embryos were fertilized in vitro according to standard protocols177 and reared in either 0.1X Marc’s Modified Ringer’s solution or 0.75X Marc’s Modified Ringer’s solution (only post stage 9). Xenopus embryos were housed at 14 °C and staged according to Nieuwkoop and Faber178.

The basal Xenobots used in this study were derived from Xenopus embryonic ectodermal explants, also known as “animal caps”, as the starting material49–52. The basal Xenobots are the most basic version of Xenobots (no sculpting or engineering, or mixing of different tissues)32. These Xenobots were constructed manually from Xenopus Laevis embryos as described previously33. Briefly, in vitro fertilized embryos were reared at 14 °C in 0.1X Marc’s Modified Ringer’s solution (MMR) until Nieuwkoop and Faber stage 9178. At stage 9, embryos were transferred into a petri dish coated with 1% agarose made in 0.75X MMR and containing 0.75X MMR solution. Using surgical forceps, the vitelline membrane was removed from the embryos, and animal cap epidermal progenitor cells were cut out of the embryos as per previously described protocols49–52. The ectodermal explants were placed on the agarose with inside surface facing up. Over a period of 2 h the ectodermal explants round up into a spherical shape. These are then incubated at 14 °C on 1% agarose in 0.75X MMR for 7 days with daily cleaning. By day, 7 the tissue differentiates and transforms into an autonomously moving synthetic epidermal entity (basal Xenobot) used in this study. These Xenobots were not further modified and were used as mature at day 7.

Total RNA was extracted from either stage 35/36 Xenopus Laevis embryos or from autonomous moving, fully mature Xenobots (7 days post-extraction of ectodermal explants from embryos). For stage 35/36 embryos, we had 3 replicates for each experimental sample, and each sample contained 10 pooled embryos. For Xenobots we had 3 replicates for each experimental sample, and each sample contained 50 pooled Xenobots—necessary to gather enough RNA from markedly small mass of Xenobots (Fig. 2A). Total RNA was extracted using Tri-reagent (MRC, Inc.) as per the manufacturer’s protocol. Total RNA quality and quantity were measured using NanoDrop spectrophotometer (Thermo Fisher Scientific).

Total RNA was sent to the University Genomic Core. RNA quality was measured with bioanalyzer, and high-quality RNA was used for library preparation with the Illunima Stranded Total RNA with Ribo-Zero Plus. Libraries were multiplexed, and single-end 75-nucleotide sequencing was performed with 30 million reads per sample on Illumina NextSeq 550 High output. Raw read files were used for initial analysis at Bioinformatics and Biostatistics Core at Joslin Diabetes Center.

The reads were trimmed for adapter “CTGTCTCTTATACACATCTCCGAGCCCACGAGAC” and polyX tails, then filtered by sequencing Phred quality (≥Q15) using fastp179. Xenopus laevis genome sequences and gene annotation were downloaded from the NCBI genome database, version 10.1. GenomeGenerate module of the STAR aligner180 was used to generate the genome indexes. STAR aligner option was set to sjdbOverhang = 75 for 76-bp reads, as was ideal. Adapter-trimmed reads were aligned to the genome using STAR aligner with the two-pass option. Reads were mapped across the genome to identify novel splice junctions in the first-pass. These new annotations were then incorporated into the reference indexes, and reads were re-aligned with this new reference in the second pass. Gene expression was estimated from the gene alignments using RSEM tool for accurate quantification of gene and isoform expression from RNA-Seq data181. Low-expressing genes were filtered out by only keeping genes that had counts per million (CPM) more than 0.6 in at least 3 samples. There were 28,009 genes after filtering. Counts were normalized by weighted trimmed mean of M-values (TMM)182. Voom transformation183 was performed to transform counts into logCPM where logCPM = log2(106∗count/(library size∗normalization factor)). PCA analysis was performed to provide an overall view of the data. Differential gene expression analysis between groups was performed using limma184.

A schematic of the method is shown in Fig. 2A. The analysis was performed on gene count data from the 3 Xenobot pools and 3 age-matched Xenopus embryo pools in MATLAB. The 28,009 genes were ranked by the values of the means of their counts across all 6 pools of Xenobots and Xenopus embryos. Genes with a count of 0 in any of the 6 pools were removed, leaving 25,276 genes. Because the count value of a given gene in a given Xenobot pool represents the mean value of all the individual Xenobots in that pool, the standard deviation of the 3 Xenobot pools gives the standard error of the means (SE) of the Xenobots, and likewise for the age-matched Xenopus embryos. The SE is related to the number of individuals per pool (n) and the standard deviation of the individuals within the pools (σ) with the equation SE = σ/sqrt(n). By multiplying the SE by sqrt(n) (where n = 50 for Xenobots, n = 10 for age-matched Xenopus embryos), σ can be calculated. Dividing σ by the mean count value of the 3 pools gives the coefficient of variation (CV) for Xenobots (CVX) and for age-matched Xenopus embryos (CVE). The CVs for all genes in the Xenobots and age-matched embryos were plotted as histograms with the area under the curves normalized to 1 (Fig. 2B) and the two distributions compared with a Wilcoxon Rank Sum test in MATLAB, returning p = 0, indicating highly significant difference between the distributions. The ranked gene list was then split into equal-size bins (percentiles), and the fraction of genes in each bin for which CVX > CVE was calculated and plotted in Fig. 2C, with black line indicating fraction of 0.5. The fraction of all genes for which CVX > CVE was found to be 0.9606 and was plotted as the red line.

To assess whether any bin values were outside the range of values expected for random shuffles of the list, a permutation test was used to determine the typical distribution of bin values around the red line. First, the order of CVX–CVE gene pairs list was shuffled, and the bin analysis was performed on the shuffled list to determine a set of bin fractions. This was repeated 1000 times for different shuffles to produce a distribution of bin fraction values for each bin. To determine statistical significance of the true bin values’ differences from the mean (red line), each true bin fraction value was compared to the distribution. The p value was defined as the proportion of the distribution that was further (in absolute value) from the distribution mean (red line) than the true bin fraction. Bins with p values of p < 0.05 were deemed statistically significant and were colored dark blue; otherwise, bins were colored light blue, revealing significant variation of bin fractions across gene count percentiles. To generate the list of top 10 most variable Xenobot genes CVX values were used. For the list, genes for which any transcriptome in the group of 3 replicates had a count value of 0, and genes for which all count values in the 3 replicates were <10, were excluded.

To understand the biological processes and pathways enriched in RNA-seq comparisons, we performed functional enrichment analysis using the Database for Annotation, Visualization and Integrated Discovery (DAVID)81,82. Gene lists of interest (significantly differentially expressed genes with false-positive genes subtracted) were uploaded to the DAVID Analysis Wizard. Functional Annotation Clustering was performed to identify significantly enriched functional annotations, and these annotations were grouped into related clusters based on shared genes. Stringency classification was set to the default recommendation of “Medium”. Each individual annotation was assigned a p value, and each cluster of annotations received an enrichment score, which is the geometric mean (in −log scale) of members’ p values in a corresponding annotation cluster. Enriched functional clusters were ranked based on their group enrichment score.

We performed a network analysis for identifying active functional biological modules83,84. We combined gene expression and interaction data, extracted the xenobot protein-protein interaction, and applied network embedding followed by clustering, similarly to85,86. We used the MNMF network embedding and clustering algorithm185 and imposed a number of clusters equal to 15. To construct the synthetic proto-organism network, we filtered the genes of interest, found the human orthologs in the different conditions using HCOP186. Then, we extracted the corresponding network using the STRING database187. The detected functional modules after network embedding and clustering are enriched using g:profiler87.

We used phylostratR to perform the phylostratographic analysis of xenobot transcripts188. The focal species was set to “8355” for Xenopus Laevis. This package automates several key processes in evolutionary analysis (1) it constructs a clade tree from species listed in UniProt, aligned with the latest NCBI tree of life; (2) it prunes the clade tree such that it retains a phylogenetically diverse selection of representatives for each phylostratum; (3) it compiles a database of protein sequences from hundreds of species from the constructed clade tree and sourced from the UniProt Proteome database (in total 329 species in our study, including the yeast and human proteomes that we added manually) (4) it applies similarity search by conducting pairwise BLAST between the proteins encoded in the focal species (here the Xenopus Laevis) and all proteins from all target species within the clade tree; (5) it determines the “best hits” and infers homology for each gene of the focal species against each target species; (6) it assigns each gene to a phylostratum correlating with the oldest clade for which there is an inferred homolog. Genes found exclusively in the focal species are identified as orphan genes and are categorized under the phylostratum “Xenopus Laevis”. The phylostratas we selected are: All living organisms (Eubacteria, bacteria and their descendants), Eukaryota, Opisthokonta, Metazoa, Eumetazoa, Bilateria, Deuterostomia, Chordata, Vertebrata, Gnathostomata, Euteleostomi, Sarcopterygii, Tetrapoda, Anura, Xenopus, Xenopus Laevis. We then mapped the age of genes with our xenobot overexpressed genes in the different conditions and the control.