Embryos assist morphogenesis of others through calcium and ATP signaling mechanisms in collective teratogen resistance

A number of questions remain open. First, we do not know how this effect operates in the wild or how the many additional factors of a complex natural environment will impact CEMA. Thus, while we have shown that embryos can cooperate in groups to resist teratogenesis, we do not yet know how much this affects evolution and adaptive fitness in nature. We also do not know the limits of this effect across types of injury. We do not claim it is universal, and many more perturbations beyond the three tested above will need to be tested to see how broadly general this effect may be. It would be interesting to see what other types of developmental perturbations can cause signal generation and propagation. Subsequent experiments will look at gene expression changes in individual animals of a cohort, paralleling the strategy of doing single-cell RNAseq in single bodies156–158, to help understand the scaling of instructive cues and their transcriptional responses.

Other avenues for future work are now open. This needs to be investigated in other systems, particularly mammals. It is already known that calcium waves exist across mammalian tissues92,94, and that mammalian cells in culture do much better in groups than they do alone159,160, but we don’t know if this extends to cohorts in embryonic mammalian development. The general field of Allee effects35–37 in novel embodiments is a fascinating area for future investigation, especially as it may connect with other instances of horizontal transfer of influence via substrates of different scales, from exosomes161–163 to whole tissues or even organs transplanted across bodies164–170. We hope that the above dataset and unconventional assay form the basis for an integrated approach to understanding the robustness of collective decision-making across scales, adding to the growing literature on collective problem-solving18–20,171–174.

There are also potential connections to evolutionary questions related to group selection175–177, which may be enriched by a better understanding of what group dynamics contribute to the developmental fitness of each member of a community. It is likely that interesting evolutionary developmental biology analyses could emerge from a broad study of CEMA, impacting basic questions about the information flows that determine embryonic outcomes178,179. In general, these types of inputs into developmental outcomes stretch the concept of epigenetics via environmental influences. Moreover, cross-embryo beneficial influences challenge the assumption that non-genetic external factors are always due to inanimate physical aspects of the environment, conspecific competition, or harmful exploitation by other biotic agents180.

Taken together, these data reveal an essential role for inter-embryo signaling during morphogenesis, adding to the growing body of work indicating that the relationship between genomic information and anatomical outcome is not as straightforward as expected. For example, planaria accumulate many mutations through somatic inheritance but are able to regrow any missing part with 100% fidelity to reach a target morphology181. The fact that organisms can reach their correct morphology despite variation in genetic material shows that there is a gap in understanding all the informative sources of morphology. It is essential to continue to elucidate diverse instructive inputs to understand the mechanisms of anatomical plasticity and robustness.

The discovery of the CEMA effect may have several practical impacts. The first concerns available data on the developmental toxicity of various agents. Our results indicate that the percentage of defects (i.e., teratogenic potential) identified in assays such as FETAX182–185 are unwittingly adjusted for CEMA—the degree of teratogenicity reported is what is seen after an unknown level of embryonic assistance has taken place (since these assays are almost never done on singletons). In other words, because most studies do not compare (or sometimes even state) the size of the cohorts, we rarely know the actual effect of a given agent—only its effect after possible CEMA mechanisms have had a chance to improve it. This suggests that it is essential to state the cohort size in teratogenicity assays, and even to re-do many of the most important studies using different embryo group sizes, to improve the transferability of the findings to human embryos where the cohort size is much smaller than in Xenopus or Zebrafish models.

The more positive implication concerns how understanding CEMA might improve biomedicine. It is tempting to speculate that, having understood how a group of agents signals each other to establish a healthier outcome, we may someday be able to artificially trigger that response in therapeutic settings. While much more work needs to be done to properly evaluate the potential of this effect for biomedicine, it is clear that beneficial, instructive cross-embryo communication, understood broadly, is an exciting phenomenon that could shed light on evolutionary fitness of developmental mechanisms and could perhaps be hacked to address urgent biomedical needs in birth defects and other disorders of morphogenesis186,187.

Animal care was done in compliance with and approval from the Institutional Animal Care and Use Committee (IACUC) under protocol number M2023-18 of Tufts University.

Xenopus laevis tadpoles were reared in 0.1× Marc’s Modified Ringers solution (MMR), pH 7.8 using standard procedures188 and were staged according to Neiuwkoop and Faber189. All embryos (pooled from separate mothers and then randomly divided between treatments) were raised at 14 °C. Feeding stage animals were fed 3 times a week with Sierra Micron powdered diet, and media changes were performed M, W, and F. Tadpoles were raised in either small dishes or large dishes depending on the cohort size.

Thioridazine (thioridazine-HCl, Sigma-Aldrich) was dissolved in deionized water at 1 mM and frozen in single-use aliquots to prevent continuous thawing. Neurula stage embryos that had been reared at 14 °C prior to exposure were treated with 90 µM thioridazine in 0.1x MMR from NF stage 12.5 to 26. During the treatment, experimental and control embryos were kept at 18 °C. Post-exposure they were returned to 14 °C until they reached scoring stage 45. A ratio of 1 mL of media to 1 embryo was kept constant across group sizes (therefore 40 embryos were raised in 40 mL and 120 embryos were kept in 120 mL).

Forskolin (Forskolin, Tocris) was dissolved in dimethyl sulfoxide at 10 mM, aliquoted into single doses, and frozen until use. Stage 10 embryos were exposed to 5 µM forskolin in 0.1x MMR. Animals were reared at 14 °C and media was refreshed three times a week (M, W, F) until animals reached scoring stage 45.

Nicotine (Sigma-Aldrich) exposure occurred from stage 11 to stage 35 at 0.1 mg/mL nicotine in 0.1X MMR, with media/drug refreshed every other day. Animals were kept at 18 °C during treatment. Post-treatment, animals were rinsed with 0.1x MMR twice, moved into fresh 0.1x MMR, and kept at 14 °C until stage 45.

Suramin (Suramin sodium salt, Sigma-Aldrich) was dissolved in dimethyl sulfoxide at 10 mM, aliquoted into single doses, and frozen until use. Exposure occurred with 100 µM suramin in 0.1x MMR conducted on neurula stage embryos (NF stage 12.5–26) that had been reared at 14 °C prior to exposure. During the treatment, experimental and control embryos were kept at 18 °C. After the exposure, they were returned to 14 °C until they reached scoring stage 45. Calcium wave experiments were performed on pre-neurula stage embryos (NF stage 10–12). Embryos were exposed to 100 µM suramin in 0.1x MMR for 30 min prior to treatment at 18 °C. Embryos were maintained in suramin during the measurement of spontaneous and post-injury calcium activity.

Tadpoles used for morphometric analysis were imaged with a Nikon SMZ1500 microscope with a Retiga 2000R camera and Q-capture imaging software. Landmarks for morphometric analysis were based on reproducibility across tadpoles. ImageJ was used to measure the width between two points (head width). For head width, landmarks were the outermost points of the head (generally near the middle of the eyes), marking the diameter and base of the branchial arches.

mRNA for a chimeric construct that is dominant negative for Kir6.1 was synthesized using standard message machine kits (Life Technologies) and stored at −80 °C until used. Embryos were transferred to 3% Ficoll solution before being microinjected. Pulled capillary needles were used with bubble pressure between 55–60 kPA and injection pressure set to 140 kPA. Injection time was set at 100 ms and at the 2-cell stage, 2 out of 2 cells were injected. Immediately after injection, embryos were moved to fresh 3% Ficoll plates and left to recover for 1 to 2 h. After that timeframe, half the media was poured out and filled with 0.1x MMR for another hour. Afterward, embryos were washed twice in 0.1x MMR and moved to a 14 °C incubator. Similarly, mRNA for the genetically encodable fluorescent calcium reporter GCAMP6S was synthesized from linearized template DNA by the message Sp6 in vitro transcription kit (Thermo Fisher) using previously described methods190,191. Roughly 2 nl of 300 ng/ml mRNA was delivered at the 4-cell stage to 4 out of 4 cells by microinjection.

GCAMP6S-injected embryos (NF stage 10–12) were loaded in groups of 2 into custom-machined acrylic holders with channels (4.6 × 1.6 × 2.5 mm). A thin layer of mineral oil (Sigma-Aldrich) was added on top of each channel to prevent drying during imaging sessions. Using a ZEISS Axio Zoom.V16 microscope and frame rate of 1 image acquired every 2 s, 10 min of spontaneous GCAMP6S activity was captured for all embryos at baseline. Localized mechanical injury was induced in 1 of the 2 embryos by a pulled glass capillary needle. Using identical imaging parameters, GCaMP6s expression was captured for all embryos from the time of injury for a duration of 20 min.

In Fiji (ImageJ), image stacks including baseline and post-injury activity were cropped to separate image stacks for each embryo (injured and receiver) and each stack was registered using the descriptor-based time series registration plug-in. A circular region of interest was selected to encompass each embryo and the mean gray value was output for that region across all frames acquired. The mean gray values after injury were normalized to the baseline activity collected for each embryo.

The peak of activity for a given region was defined as the maximum normalized signal at any time during post-injury data capture. The distance between any 2 regions of interest was calculated by finding the magnitude of the vector from the injury site to the closest point on the receiver embryo. The time of calcium wave transfer between embryos was taken as the time from injury induction to the time that the mean calcium signal surpassed the maximum of spontaneous activity within an embryo. To avoid detecting aberrant spikes, a transfer between embryos was only considered complete if the mean calcium signal was maintained above spontaneous levels for longer than 1 min. The cell-to-cell speed within a single embryo (injured embryo & receiver embryo speeds) was measured by selecting a cell at the start of a calcium wave and a cell at the end of the calcium wave. The time from cell-to-cell was calculated by measuring the peak-to-peak time delay between each cell. Speed for intra- and inter-embryo measurements was defined as the time of calcium wave transfer divided by inter-region distance.

All models were built using the MATLAB software (9.12.0.1884302, R2022a). Each ECA was created with 149 cells with a full runtime of Ncells/2, a common configuration82,83. To simulate development, ECAs were given the GKL rule to solve the majority problem, with an initial configuration of 60% 1 s and 40% 0 s84. An embryo was considered successfully developed with no defects if all cells had converged to 1 before the end of simulation time. To simulate embryo groups, they were developed in a square configuration (2 × 2, 3 × 3, 10 × 10, for example) from values ranging from singletons (1 embryo) to large groups (18 × 18, or 324 total embryos). Tested values were chosen to emulate the experimental data.

To simulate health and noise (teratogen), each embryo was given an initial health value of 1. At each timestep there is a chance that noise, pnoise, could affect any embryo in the group, altering its health by a percentage, ndec, parameterized to 0.7. Then, each embryo affected by the noise will update its health value based on a weighted average of its neighbors and neighbor’s neighbors. The embryo’s weight is 1, its neighbor’s weight is 1, and its neighbor’s neighbors’ weight is 0.25. These values were chosen to simulate spatially close mixing as well as diffuse, slightly distant signaling. This signaling was vital to the CEMA effect, as when the weights were turned low (<0.01) the collective never recovered unless noise was turned down to a substantial degree (pnoise = 0.6) (see Supplementary Fig. 4A). After, each embryo affected by noise will again update its health based on supportive signaling, sampled from the previous timestep. This health update is also calculated based on a weighted average, with the weights of the embryo’s own signal proportional to (1-health), as is the signal from its neighbors. The weights of the neighbor’s neighbors are weighted as 0.25*(1-health), again to simulate diffuse signaling. This weighted average is then applied to the embryo’s health value following the formula:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{{\mathrm{health}}}}}}}\,{{{{{{\mathrm{value}}}}}}} +={{{{{{\mathrm{current}}}}}}}\,{{{{{{\mathrm{health}}}}}}}\,{{{{{{\mathrm{value}}}}}}}*{{{{{{\mathrm{weighted}}}}}}}\,{{{{{{\mathrm{support}}}}}}}\,{{{{{{\mathrm{signal}}}}}}}$$\end{document}healthvalue+=currenthealthvalue*weightedsupportsignal

After this, an embryo then evolves according to the GKL rule if health is above 0.5, or the random update rule if health is equal to, or less than, 0.5. This threshold was chosen as the intermediate value between 0 and 1, however, other values were tested (see Supplementary Fig. 4B). At higher values (0.75) embryos deteriorated too quickly for recovery, and for lower values (0.25), recovery became trivial, even for smaller conspecific groups. For parameterization, noise values from 0.0 to 1 were tested in steps of 0.1, and we found that noise of 0.8 best reflected the data found experimentally. For ndec, parameterization found that 0.7 best fit the experimental data. Each simulation was run 20 times to create a sample space, and 99% confidence intervals were calculated for each noise value across a number of embryos.

To simulate the experimental design when embryos not affected by the teratogen were introduced into a group that had previously been exposed, we artificially locked half of the population to 1, meaning they stayed perfectly healthy. With this, they were not allowed to participate in the health signaling.

At stage 35, embryos were sacrificed for an rRNA depletion study. A sample consisted of 15 pooled tadpoles, and each was repeated 3 times. Tissue was extracted using TRIzol (Thermo Fisher Scientific) as per the manufacturer’s protocol, and total RNA quality and quantity were assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific).

RNA was sent to the Tufts Genomic Core where RNA quality was assessed via bioanalyzer, and high-quality RNA was used for library preparation with the Illumina Stranded Total RNA with Ribo-Zero Plus. Libraries were then multiplexed, and an rRNA depletion run using single-end, 75-nucleotide sequencing was performed on Illumina HiSeq 2500. Raw read files were sent to the Bioinformatics and Biostatistics Core at Joslin Diabetes Center.

Reads were trimmed for adapter “CTGTCTCTTATACACATCTCCGAGCCCACGAGAC” and polyX tails, then filtered by sequencing Phred quality (>=Q15) using fastp192. Adapter-trimmed reads were aligned to the Xenopus laevis rRNA genomic sequence (version 10.1) from the NCBI nucleotide database using Bowtie2193 and unmapped reads were removed using samtools194. Adapter-trimmed reads were aligned to the genome using a STAR aligner with the two-pass option. The RSEM tool was used to estimate the gene expression from the gene alignments. Low-expressing genes (threshold of 1.8 counts per million in at least 3 samples) were filtered out, leaving a total of 21,720 genes after filtering. Counts were then normalized by the weighted trimmed mean of M-values (TMM)195. The counts were Voom transformed196 into logCPM (logCPM=log2(106∗count/(library size∗normalization factor))). Two surrogate variables (SVs) were identified and constructed197. Batch and the SV effects were adjusted and a principal component analysis (PCA) was performed. Differential expression analysis between groups was performed using limma198. Pathway analysis was performed by testing the over-representation of the differentially expressed genes (FDR < 0.1) in the Gene ontology (GO) terms/pathways using R package clusterProfiler199. Enriched pathways were identified at a significance threshold of adjusted p-value (FDR) < 0.1. In the gene ontology dot plots, the size of the dots reflects the gene ratios (number of significant genes associated with the GO term/total number of significant genes associated with any GO term), and the adjusted p-value (FDR) reflects the significance.

All statistical analyses were performed using Prism 9. To achieve statistical power, biological replicates (N) were conducted 3–6 times with n > 50 embryos for each treatment unless otherwise noted. Data across various iterations were pooled and analyzed by non-parametric t-test (for two groups) or ANOVA (for more than two groups).

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Supplementary Information Peer Review File Description of Additional Supplementary Files Supplementary Movie 1 Supplementary Movie 2 Supplementary Movie 3 Supplementary Movie 4 Reporting Summary