The Morphological, Behavioral, and Transcriptomic Life Cycle of Anthrobots

Abstract

Fascinating aspects of morphogenetic and behavioral plasticity of living material are revealed by novel constructs that self-construct from genetically wild-type cells. Anthrobots arise from cultured adult human airway epithelial cells, developing, becoming self-motile, and acquiring neural repair capabilities without exogenous genetic circuits or inorganic scaffolds. Progress in bioengineering and regenerative medicine depends on developing a predictive understanding of collective cell behavior in novel circumstances. Toward that end, here a number of life cycle properties of Anthrobots, including their morphogenesis, maturation, and demise, are quantitatively characterized. A self-healing capacity and a remarkable reduction of epigenetic age upon morphogenesis are uncovered. Transcriptomic analysis reveals that assembling into Anthrobots drives a massive remodeling of gene expression relative to their cellular source, including several embryonic patterning genes, and a shift toward more evolutionarily ancient gene expression. These data reveal new aspects of engineered multicellular configurations, in which wild-type adult human cells self-assemble into an active living construct with its own distinct transcriptome, morphogenesis, and life history.

Keywords: biobots; repair; synthetic morphology; transcriptomics.

Figure 1

Anthrobots have different modes of self‐construction that may account for the emergent morphotypes. A) Principal components analysis (PCA) cloud consisting of 28 Anthrobots’ growth time course biomass data, evaluated with both morphological and temporal indices. The three orthogonal clusters identified by the unsupervised clustering algorithm are labeled as groups 1, 2, and 3, indicating three major fates respectively: dormant, merger expander, and monoclonal expander. The indices forming the PCA dimensions are as follows. The “merge_events” index is the master metric that distinguishes between Anthrobot formation via monoclonal expansion versus cluster merging: if a spheroid is undergoing monoclonal expansion, this number would be 1, while if the final spheroid is a result of multiple spheroids merging together, this number would denote the number of times a merge event has occurred. The “fin_area” measures the area of the final spheroid, while “sd_area” quantifies the average standard deviation of area across entire growth process. B) These three modes of growth are further described with the associated boxplots (n = 9 for group 2, n = 3 for group 3, and n = 16 for group 1). C) Representative examples from the PCA, as labeled on the cluster plot. c1) An example dormant bot. Notice volume stays the same. Scale bar 10 µm. c2) An example of hybrid Anthrobot formation through a merge event. Notice the heterogeneous structure in comparison to c1 and c3. Scale bar 10 µm. c3) An example monoclonally expanding, imaged bidaily. Scale bar 50 µm. p‐value range of 0–0.0001 corresponded to ****, 0.0001–0.001 corresponded to ***, 0.001–0.01 corresponded to **, 0.01–0.05 corresponded to *, and 0.05–1 corresponded to ns. All significance tests were evaluated at an alpha value of 0.05. Unless otherwise specified, the alternative hypothesis was always two‐sided for Wilcox tests.

Figure 2

Early stage Anthrobots display onset markers of embryonic development. As NBHEs progress from progenitor cells to Anthrobots, we observe gene patterning characteristics of mammalian germ layer development and axis formation. A) Example morphologies of “progenitor cell,” at the beginning of Anthrobot growth time course, and “day 0 bot,” meaning it is on the last day of bot formation time course and day 0 of its spheroid (i.e., bot state) time course. B) PCA of clustering of Anthrobots across different life stages. Progenitor and day 0 bot stages are shown with dashed circles. C) The difference in gene expression between the progenitor versus day 0 stage, showing significant differences in the transcriptome, with notable genes explained further in panel (D) (p < 0.005; n = ≈1 million progenitor cells). D) Functional genes for germ layer formation (FOXA1, FOXA2, FOXC1, BMP7, and FOXP2) as well as axis formation (BMP7, HESX, and SHH) were observed to display expression profiles of early embryonic development. E) Histogram of the phylostratigraphic analysis of the Anthrobots and progenitor cells in the different conditions, showing the number of genes expressed with log counts per million (CPM) > 1 in the different conditions. F) Number of overexpressed DEGs with logFC > 2 in the different conditions. For (E) and (F), the X‐axis shows the evolutionary ages of ancient genes (here “All living organisms” and “Eukaryota,” the former corresponds to Eubacteria, bacteria, and their descendants). The Y‐axis shows the gene counts. See Supplement S1 in the Supporting Information for additional details. N = ≈1800 bots for day 0, 1000 bots for day 10, and 600 bots for day 25. p‐value range of 0–0.0001 corresponded to ****, 0.0001–0.001 corresponded to ***, 0.001–0.01 corresponded to **, 0.01–0.05 corresponded to *, and 0.05–1 corresponded to ns. All significance tests were evaluated at an alpha value of 0.05. Unless otherwise specified, the alternative hypothesis was always two‐sided for Wilcox tests.

Figure 3

Anthrobots display different behaviors during the eversion process. A) Sample Anthrobot eversion time course. Day 0 Anthrobot going through polarity reversal via eversion across 48 h. B) PCA cloud consisting of n = 27 Anthrobots’ eversion time course biomass data, evaluated with both morphological and behavioral indices. The two orthogonal clusters identified by the unsupervised clustering algorithm are labeled as groups 1 and 2, indicating two major modalities of eversion. Plotting high‐dimensional cloud 2D. C) These two modes of eversion are further described with the associated boxplots (n = 12 for group 1, n = 15 for group 2). More specifically, while a majority of Anthrobots stay in place during eversion (group 2), which would be the expected behavior, a smaller subset of bots seems to be displacing during polarity reversal (group 1) as shown on graph i (p = 0.00027). This motility is likely caused by disorganized forces generated by radical anatomical activity, causing large variability in the heading of the displacing bots as shown in graph ii (p = 0.00018). The fact that there is no difference in terms of the metrics that measure organized behavior such as gyration (graph iii) or straightness (graph iv) between these two modes of eversion further supports this explanation. Straightness p = 0.051, gyration p = 0.8. D) Representative examples from each group across time. The individual bots are marked on the PCA cloud in panel (B) as botD1 and botD2. E) Collection behavior demonstrated by one of the motile eversion bots (group 2). The residual filament in cell culture media (shown with an arrow) is collected by a group 2 bot that ultimately encases the filament debris during eversion. Asterisks indicate statistical difference to p < 0.05. p‐value range of 0–0.0001 corresponded to ****, 0.0001–0.001 corresponded to ***, 0.001–0.01 corresponded to **, 0.01–0.05 corresponded to *, and 0.05–1 corresponded to ns. All significance tests were evaluated at an alpha value of 0.05. Unless otherwise specified, the alternative hypothesis was always two‐sided for Wilcox tests.

Figure 4

Late‐stage Anthrobots display markers of maturing embryonic development. As NBHEs progress from apical‐in state on day 0 to apical‐out state on day 10, we continue to observe gene patterning characteristics of mammalian germ layer development and axis formation. Colors green, blue, red, and yellow, respectively, represent basal cells, nuclei, tight junctions, and cilia. A) Example morphologies of “day 0 bot,” at the beginning of Anthrobot eversion time course, and “day 10 bot,” meaning it is at a mature motile state (n ≈1800 bots for day 0, n ≈1000 bots for day 10). B) PCA of clustering of Anthrobots’ transcriptome across different life stages. Day 0 and day 10 bots are shown with dashed circles. C) The difference in gene expression between the day 0 versus day 10 stage, showing significant differences in the transcriptome, with notable genes explained further in panel (D). D) FOXA1 and FOXA2, key transcription factors for lung epithelial cell proliferation and differentiation in endoderm, show no difference of expression in day 0 versus day 10 bots. Similarly, genes like FOXC1 and BMP7 that are associated with mesoderm formation show decrease in expression. Conversely, an increase in the gene FOXP2, ectoderm marker, is seen as upregulated. Furthermore, for axis formation, while the dorsal–ventral patterning marker BMP7 gets downregulated in this later developmental stage, the anterior–posterior marker HESX stays steady and the previously silent left–right patterning marker SHH becomes upregulated on day 10. p‐value range of 0–0.0001 corresponded to ****, 0.0001–0.001 corresponded to ***, 0.001–0.01 corresponded to **, 0.01–0.05 corresponded to *, and 0.05–1 corresponded to ns. All significance tests were evaluated at an alpha value of 0.05. Unless otherwise specified, the alternative hypothesis was always two‐sided for Wilcox tests.

Figure 5

Anthrobots display robust behavioral longevity and reversal of epigenetic age. A) Weekly counts of baseline behavioral activity (defined as any cilia‐induced movement, regardless of the amount displaced) stay steady despite the shrinking population due to aging. B) Toward the end of an Anthrobot lifecycle, between weeks 4–8, while there is visible degradation of bot structural coherence (see Figure 6), the persistence of baseline behavior stays at a steady frequency. Population‐level baseline behavior outlives individual bots. C) Cells harvested from a 21 years old donor were differentiated into Anthrobots, which were then collected at day 10 and day 25 time points (n ≈ 1000 bots for day 10, and n ≈ 600 bots for day 25 time points). Methylation clock studies revealed passage 1 primary cells’ tissue age as 23.7 (mean value) years old (n = 12). Both against this baseline and against the recorded age of the donor, differentiation into Anthrobots introduced significant methylation clock reversal, with day 10 bots being read at the mean value of 19.1 years old (n = 6, p = 0.0015), while the day 25 bots being read at the mean age of 20.0 (n = 3, p = 0.0087). Contrarily, an air–liquid interface differentiation of the same cells into 2D airway tissue (ALI tissue) increased tissue age to a mean of 25.7 years old (n = 3, p = 0.0025 versus day 10 Anthrobot time point). D) Day 4 bot's healing process from a hypodermic needle injury (see Supplement S3 in the Supporting Information for additional details).

Figure 6

Regardless of initial size, Anthrobots degrade at a steady rate, resulting in large bots living longer. A) PCA analysis showing 36 bots’ cluster analysis based on morphological indices characterizing the change in size across a 30 days timeframe. B) Bots that start with a i) large initial area (group 2, n = 10) compared to bots that start with a smaller initial area (group 1, n = 26) show the same ii) relative percent change across their degradation time course, maintaining the size difference iii) at the end of the degradation process. Regardless of initial morphology, all Anthrobots are subject to a similar rate of aging, resulting in large bots living longer than smaller bots. Panels (i) and (iii) use arbitrary units in terms of pixel square. C) Bot C1 from group 1 (marked on the PCA) has a smaller initial area compared to bot C2 from group 2 with a larger initial area. While bot 1 completely degrades within the same time frame, bot 2 still stays intact, suggesting a positive correlation between bots’ initial morphology and their life span as opposed to a finite time span across all population. Scalebar 50 µm. Asterisks indicate statistical significance to p < 0.05. p‐value range of 0–0.0001 corresponded to ****, 0.0001–0.001 corresponded to ***, 0.001–0.01 corresponded to **, 0.01–0.05 corresponded to *, and 0.05–1 corresponded to ns. All significance tests were evaluated at an alpha value of 0.05. Unless otherwise specified, the alternative hypothesis was always two‐sided for Wilcox tests.