The Morphological, Behavioral, and Transcriptomic Life Cycle of Anthrobots
By characterizing new ways in which adult somatic cells can reboot their multicellularity, and understanding the way these wild‐type cells navigate transcriptional, morphological, and behavioral spaces, novel aspects of the cell–body relationship are revealed. Novel transcriptomic profiles, driven by configurations and lifestyle, as well as behaviors not exhibited in the standard species‐specific morphological trajectory, have implications for basic evolutionary, cell, and developmental biology and the appearance of novelty and evolvability.[ 49 ]
The NHBE cells were seeded into a Matrigel suspension according to the previously established protocol.[ 3 ] Once in the Matrigel suspension, select NHBEs (15–20 beacons/points of interest) were imaged using an EVOS M7000 (AMF7000) once every 3–4 days using a 20× objective in both Brightfield and Red Fluorescent Protein (RFP) until they had been 14 days in the Matrigel, at which point they were removed from the Matrigel suspension and used in downstream applications.
RNA was extracted from Anthrobots and NHBE cells under different time points and experimental conditions. The extractions were performed using Invitrogen's TRIzol Reagent user guide as a reference. For the NHBEs, the only additional step was to trypsinize the cell layer with. 05% Trypsin for 3–4 min until the cells had been suspended in their media. Once the material (Anthrobots and NHBEs) was in suspension, it was collected into 15 mL conical tubes and centrifuged at 300 g for 3 min. The supernatant was then removed, and the material was resuspended in 1 mL of TRIzol and pipetted up and down to homogenize. Next, 0.2 mL of chloroform was added and then the sample was centrifuged for 15 min at 12 000 g and 4 °C. After centrifugation, the aqueous phase containing the RNA was transferred to a new tube. 0.5 mL of isopropanol was added to the aqueous phase and incubated at room temperature for 10 min. After incubation, the sample was centrifuged for 10 min at 12 000 g and 4 °C. The RNA formed a small, white pellet on the bottom of the tube. Next, the supernatant was discarded, being careful not to disturb the RNA pellet. The pellet was then resuspended in 75% ethanol, vortexed briefly, and centrifuged again for 5 min at 7500 g and 4 °C. The supernatant was then discarded, and the pellet was left to dry for between 5 and 10 min. After drying, 50 µL of RNAse‐free water was added to the pellet, and the sample was incubated on a heat block at 55 °C for 10 min. Finally, the RNA concentration in the sample was measured using a NanoDrop One from Thermo Fisher Scientific.
RNA quality was assessed via bioanalyzer at the Tufts Genomic Core, and Illumina Stranded Total RNA with Ribo‐Zero Plus was used for library preparation with high‐quality RNA. Upon library multiplexing, on Illumina HiSeq 2500, an rRNA depletion run using reversely stranded, single‐end sequencing was performed.
Reads were trimmed for adapter “GATCGGAAGAGCACACGTCTGAACTCCAGTCAC” and polyX tails, then filtered by sequencing Phred quality (≥Q15) using fastp.[ 50 ] A count table was generated by aligning reads to the human transcriptome (Ensembl version 110) using kallisto[ 51 ] and converting transcript counts to gene counts using tximport.[ 52 ] To filter out low expressing genes, genes that had counts per million (CPM) more than 0.29 in at least 3 samples were kept. This CPM threshold was used so that genes that were expected to have at least 10 counts were kept, which was a rule‐of‐thumb value. There were 22 518 genes after filtering. Counts were then normalized by weighted trimmed mean of M‐values (TMM).[ 53 ] If no normalization was needed, all the normalization factors would be 1. Here, the normalization factors were between 0.84 and 1.26. To use linear models in the following analysis, Voom transformation[ 54 ] was performed to transform counts into logCPM, where logCPM = log2(106 × count/(library size × normalization factor)). Voom transformation estimated the mean–variance relationship and used it to compute appropriate observation‐level weights so that more read depth gave more weights. To get an overall view of the similarity and/or difference of the samples, PCA was performed.
In order to achieve the phylostratigraphic analysis of the DEGs of the anthrobots in the different conditions, the evolutionary ages of all of the 19 660 protein‐coding human genes were first extracted from Litman et al.[ 55 ] Litman et al. assigned these ages to 1 of the 19 major phylostrata defined by Domazet‐Loso and Tautz.[ 56 ] These phylostrata were the following: all living organisms, (Eubacteria, bacteria, and their descendants), Eukaryota, Opisthokonta, Holozoa, Metazoa, Eumetazoa, Bilateria, Deuterostomia, Chordata, Olfactores, Craniata, Euteleostomi, Tetrapoda, Amniota, Mammalia, Eutheria, Boreoeutheria, Euarchontoglires, and Primates. The first two were selected to assess how much the anthrobots were behaving like unicellular organisms in terms of genetic expression and the number of genes expressed in the different conditions was counted.
The list of morula genes was extracted from Yan et al.[ 57 ] To do so, the genetic expression of each cell was averaged at the morula stage and only the genes with Reads Per Kilobase per Million mapped reads (RPKM) > 2 were kept. The overlap between the DEGs in the different conditions was computed with this list of morula genes. Then, the overlapped list of genes was enriched using g:Profiler[ 58 ] (see also Supplement S4 in the Supporting Information).
The imaging for the polarity reversal was done using Anthrobots recently dissolved from their Matrigel suspension in order to capture their inversion from apical‐in to apical‐out. The setup used involved a Bioscience Tools Stage Incubation system consisting of a CO2 Mixer (CO2‐MI), a top and bottom stage incubator (TC‐MWPHB), a temperature controller (TC‐1‐100i), and a humidifier (CO2‐500 mL). This was set up on an EVOS FL Auto microscope (AMAFD1000) used for all the imaging. After the imaging setup, the Anthrobots were removed from their Matrigel suspension following the protocol in the original Anthrobots paper (Motile Living Biobots Self‐Construct from Adult Human Somatic Progenitor Seed Cells). 200 µL of bronchial epithelial differentiation media (BEDM) was added to the inner wells of a 96‐well ultralow attachment round‐bottom plate from Corning (Cat no. 7007). The outer wells were filled with sterile DiH2O, which were found critical to prevent evaporation in the stage incubation setup. Once the plate was filled with media, the bots were moved to individual wells in a 96‐well ultralow attachment round‐bottom plate from Corning (Cat no. 7007). With one bot per well, the plate was placed into the stage incubation setup. The wells were imaged at 10× every hour for 72 h.
To find if any patterns in the behavior of bots during the polarity reversal stage could be computationally identified, Fiji was used to track the change in position of the bot for each frame using the AnalyzeParticles tool. After getting a coordinate analysis of each of the bots, a few variables were generated from the data. Initially, to determine the trajectory of the bot, the linear distance, speed, heading, and angular speed were estimated at each position to predict the coordinates of the following position. Of these variables, the mean, median, and standard deviation of the angular speed, linear distance, and linear speed were taken over the entire trajectory of the bot to try to ascertain its behavior over its entire trajectory. These 14 variables were part of the variables used to make the PCA. Two metrics were then derived to describe each trajectory: i) a “straightness” index computed as 1 minus the circular variance of the headings during the block (a value of 1 indicated a perfectly straight line) and ii) a “gyration” index computed as 1 minus the circular variance of the angular speed during the block divided by the circular variance of the same angular speeds and their additive inverse, which helped in taking into account the magnitude of the angular speeds themselves (a value of 1 indicated a trajectory following a perfect circle). These 14 variables were then used in a PCA, which had ≈63% variation in the first two components. Hierarchical clustering was then employed with the WardD2 method to derive two clusters.
DNA was extracted from Anthrobots and NHBE cells using a QIAGEN DNeasy Blood and Tissue Kit with its protocol as a reference. All buffers and mixtures were provided by the kit, saved for phosphate‐buffered saline (PBS). For the NHBEs, the only additional step was to trypsinize the cell layer with 0.05% Trypsin for 3–4 min until the cells had been suspended in their media. Once the material (Anthrobots and NHBEs) was in suspension, it was collected into 15 mL conical tubes and centrifuged at 300 g for 3 min. The supernatant was then removed, and the material was resuspended in 200 µL of PBS. Then, 20 µL of Proteinase K was added to the material in suspension. Next, 200 µL of AL Buffer was added to the suspension and was mixed by vortexing. Next, 200 µL of 95–100% ethanol was added to the suspension, and the solution was again mixed by vortexing. The mixture was then pipetted into a spin column placed inside a 2 mL collection tube and centrifuged at 6000 g for 1 min. Then, after discarding the flow‐through, 500 µL of Buffer AW1 was added, and the column was centrifuged again at 6000 g for 1 min. Then, the flow‐through was discarded, and 500 µL of Buffer AW2 was added to the column. The column was then centrifuged for 3 min at 20 000 g; this centrifugation was repeated to dry the column. After the dry spin, the column was placed in a fresh 1.7 mL tube, and 200 µL of AE Buffer was added. It was let to sit for 1 min before centrifuging for a final time at 6000 g for 1 min. Afterward, the yield of the DNA collected in the 1.7 mL tube was measured using a NanoDrop One from Thermo Fisher and stored at 4 °C until needed for downstream applications.
The DNA methylation array assays were performed using the following individual components. First was bisulfite conversion of extracted DNA using a Zymo EZ DNA methylation kit for downstream processing by either DNA methylation array technology. This was followed by MethylationEPIC BeadChip array kits, which included all the reagents required for DNA amplification, labeling, and array. The third step was to scan completed Beadchip arrays for the final DNA methylation readout. Finally, there was statistical analysis of measures of epigenetic aging for the respective species tested, an analysis of treatment versus control conditions, and an analysis of any available metadata for associations with changes to epigenetic aging. See Supplement S5 in the Supporting Information for additional details.
Anthrobots were grown according to the previously established protocol[ 1 ] and were transferred between days 2 and 11 into BEDM containing wells of a 96‐well ultralow attachment round‐bottom plate from Corning (Cat no 7007). The outer wells were filled with sterile DiH2O, which were found critical to prevent evaporation in the stage incubation setup. Next, bots were damaged with a hypodermic needle and upon withdrawing the needle, they were imaged immediately every 2 s for 10–20 min at 10× magnification using a Zeiss Axio Observer, equipped with an Axiocam 506 mono camera. The following longer term time‐lapses were acquired within a Bioscience Tools Stage Incubation system, consisting of a CO2 Mixer (CO2‐MI), a top and bottom stage incubator (TC‐MWPHB), a temperature controller (TC‐1‐100i), and a humidifier (CO2‐500 mL), using an EVOS FL Auto microscope (AMAFD1000). The wells were imaged at 10× every 30 min for between 48 and 72 h.
The Anthrobots were imaged for their degradation in the following manner. First, the Anthrobots were dissolved from their Matrigel suspension according to the previously established protocol,[ 3 ] then each well of a 96‐well ultralow attachment round‐bottom plate from Corning (Cat no. 7007) was filled with 200 µL of BEDM with Retinoic acid. After adding the media, Anthrobots were added individually to each well, resulting in 1 bot per well. The bots were imaged every 3–4 days at 20× in Brightfield using a Zeiss Axio Observer.Z1, equipped with an Axiocam 506 mono camera until they disintegrated to the point where no spheroid structure was visible. Media changes were performed every 4 days by aspirating 100 µL of media from each well and then adding 100 µL of fresh BEDM with Retinoic acid.
To computationally identify any patterns in the assembling of bots, Fiji was used to track the change in the size of the “bot area” for each frame using the AnalyzeParticles tool. After getting a CSV of each of the 28 bots, it was decided that experimental observations had three general categories needed to distinguish between 1) bots that senesced, i.e., died and stopped growing in size, 2) bots that grew completely from one cell to a bot and thus had no merges, and 3) ones that merged multiple cellular entities to form the final bot that grew to the final size. To distinguish these relationships, three variables were generated from the data, namely: i) merge_events: the number of times distinct spheroids merged to distinguish the bots that formed via mechanism (3) from those that formed via mechanisms (1) and (2); ii) fin_area: the final area of the bot at the end of the growing process which helped distinguish bots formed via mechanism (1) from those formed via mechanisms (2) and (3) since senesced bots ought to have significantly smaller area than other categories which grew to completion; iii) sdarea: the standard deviation of consecutive differences in area, i.e., the variation in the rate of change of area over time, which helped separate bots formed via mechanism (1) from those formed via mechanisms (2) and (3) since senesced bots ought to have smaller rate of change after it senesces (close to zero). This redundancy reinforced the results of the fin_area variable. These three variables were then used in a PCA, which had ≈90% variation in the first two components. Hierarchical clustering was then used with the WardD2 method to derive three clusters, which seemed to correlate in identity to the three types of bots to be distinguished.
In order to computationally identify any patterns in the bot degradation stage, Fiji was used to track the change in position, area, perimeter, and other physical parameters of the bot for each frame using the AnalyzeParticles tool. After getting a coordinate index of 36 bots, a few variables were generated from the data, focused mainly on deriving insights from the area of the main bot, so that change of area could tell about rates of degradation. These included initial area (area of bot in 1st frame), fin area (area of bot in final trackable frame), change in area (fin area − initial area), percent change (change in area/initial area × 100), SD of difference in areas (SD of change in area between successive frames), mean of difference in areas (mean of change in area between successive frames), time to disintegrate (frames till the bot is untrackable), and changerate (the maximum difference in difference of areas, akin to second derivative for area). These eight variables were then used in a PCA which had ≈82% of variation in the first two components. Hierarchical clustering was then used with the WardD2 method.
To discover the differential genes, limma was used, a R package that powers differential expression analyses.[ 59 ] Linear modeling was performed and t‐tests were moderated to detect genes that were differentially expressed between groups. The GO database has hierarchical structures so that the GO terms can be grouped by their main (parent) categories.[ 20 ] Gene Set Enrichment Analysis (GSEA) was performed in the GO database.[ 60 ] Enriched GO terms were identified at a significance threshold of adjusted p‐value < 0.05.[ 61 ] A network of the significant GO terms was constructed and the weights of edges were refined based on the Jaccard similarity.[ 62 ] Pathway communities (i.e., clusters) were then identified by multilevel modularity optimization, i.e., Louvain algorithm.[ 63 ] One GO term that had the smallest p‐value in GSEA was selected for each community. GSEA was performed in the Reactome database.[ 64 ] Enriched pathways were identified at a significance threshold of adjusted p‐value < 0.05.[ 61 ] A network of the significant Reactome pathways was constructed and the weights of edges were refined based on the Jaccard similarity.[ 62 ] Pathway communities (i.e., clusters) were then identified by multilevel modularity optimization, i.e., Louvain algorithm.[ 63 ] One Reactome pathway that had the smallest p‐value in GSEA was selected for each community.
GSEA was performed in the Gene Transcription Regulation Database (GTRD) database.[ 65 ] Enriched transcription factors were identified at a significance threshold of adjusted p‐value < 0.05.[ 61 ] A network of the significant transcription factors was constructed and the weights of edges were refined based on the Jaccard similarity.[ 62 ] Transcription factor communities (i.e., clusters) were then identified by multilevel modularity optimization, i.e., Louvain algorithm.[ 63 ] One transcription factor that had the smallest p‐value in GSEA was selected for each community.
This work was supported by Astonishing Labs, which provides a sponsored research agreement to Tufts University and operates in the regenerative medicine space.
The authors thank Bobby Brooke of the Epigenetic Clock Development Foundation for assistance with the methylation clock age characterization and analysis. The authors also thank Axel De Baat, Navneet Jawanda, and Elena Sergeeva for helpful comments on a draft of the paper; and Julia Poirier for assistance with the manuscript. RNA‐seq was performed in the Tufts University Core Facility Genomic Core, which received funding support from the NIH (Grant No. S10 OD 032203) for its Illumina NovaSeq sequencer. The authors gratefully acknowledge support via a sponsored research agreement from Astonishing Labs, and via Grant No. 62212 from the John Templeton Foundation.
Gumuskaya G., Davey N., Srivastava P., Bender A., Pio‐Lopez L., Hazel D., Levin M., The Morphological, Behavioral, and Transcriptomic Life Cycle of Anthrobots. Adv. Sci. 2025, 12, 2409330. 10.1002/advs.202409330
The data that support the findings of this study are available from the corresponding author upon reasonable request.
The data that support the findings of this study are available from the corresponding author upon reasonable request.