Mitochondrial dynamics govern whole-body regeneration through stem cell pluripotency and mitonuclear balance

Mitochondria dynamics in planarian regeneration and various cell populations

Previous studies have used vital dye combination, including SiR-DNA for DNA content, CellTracker Green (CTG) for cell volume, and Mito Tracker/PK Mito for MMP, to enrich Mitolow cells from planarian neoblasts34,35,36. To comprehensively compare the mitochondrial dynamics in various cell populations, we employed the same strategy to distinguish five cell populations based on cell cycle and mitochondrial states from wild-type planarians. Based on DNA contents, 2N and 4N cells were grouped (Fig. 1A). 4N cells were divided into two groups based on their CTG staining intensity as 4N CTGlow and 4N CTGhigh (Fig. 1B)37. Within 4N CTGhigh cells, 4N CTGhighMitolow (Mitolow) and 4N CTGhighMitohigh (Mitohigh) cells were isolated based on the level of MMP (Fig. 1C). To facilitate reproducibility, we displayed a summary of sequential plots illustrating the gating strategy for Mitolow and Mitohigh cells (Fig. S1A). To validate the Mitolow and Mitohigh cells from planarian neoblasts, we performed FACS analysis by comparing the wild type with irradiated planarians (100 Gy). The results showed that Mitolow and Mitohigh cells from SiRNeoblasts are irradiation-sensitive (Fig. S1B). We also compared the cell cycle distribution of Mitolow and Mitohigh cells with that of the nucleated cells sorted after Hoechst 33342 staining. Mitolow, Mitohigh, and nucleated cells were immediately fixed after sorting and stained with DAPI for cell cycle analysis. Both Mitolow and Mitohigh cells displayed enrichment in the S/G2/M phase (Fig. S1C). We further analyzed the ratio of Mitolow and Mitohigh cells in both X1 cells and SiRNeoblasts and found that X1 and SiRNeoblasts have similar proportions of Mitolow and Mitohigh cells (Fig. S1D). These results suggest that Mitolow and Mitohigh cells are similarly proliferating cells at S/G2/M cell cycle stages.

Fig. 1: Mitochondrial dynamics during planarian regeneration and in different cell types.
figure 1

A Density plot of FACS gating based on DNA content to isolate SiR-DNA 4N (red circle) and SiR-DNA 2N (yellow circle) cells. B Pseudocolor dot plot shows the CTGhigh (red circle, SiRNeoblasts) and CTGlow (left of the blue dot line) cell populations among SiR-DNA 4N cells. C Pseudocolor dot plot of the gating strategy for two cell populations from CTGhigh with different PK Mito Red intensity, named Mitolow and Mitohigh. D Peaks in piwi-1 intensity for the enrichment of piwi-1high, piwi-1low, and piwi-1 in each indicated cell population as determined via FACS analysis. E Colorimetric WISH images show the piwi-1 probe staining in 60 Gy-irradiated hosts at 7 days post-transplantation with Mitolow or Mitohigh cells. The numbers indicate the ratio of hosts harboring piwi-1+ neoblasts. Scale bar, 500 μm. F Rescue efficiency of Mitolow or Mitohigh cells transplanted into lethally irradiated hosts. The upper right images show a rescued animal at the indicated time points. n = 10 hosts in each of three independent experiments. Scale bar, 500 μm. G Representative Hessian-SIM images of mitochondria in Mitolow (n = 33), Mitohigh (n = 32), CTGlow (n = 39), and 2N (n = 48) cells. Scale bars, 1 μm. HK Dot plots of mitochondrial size (H), total length (I), branch length (J), and length/branch ratio (K). The quantification is shown as mean ± SEM. The unpaired two-tailed Student’s t-test was used to determine the significance of differences between two conditions. L Pseudocolor dot plots show Mitolow and Mitohigh cells with their percentages from regenerating planarians, including 0 hpa, 1 hpa, 6 hpa, 12 hpa, 1 dpa, 2 dpa, 3 dpa, 5 dpa, and 7 dpa. The quantification is shown as mean ± SEM. M Line plot shows the changes in the ratio of Mitolow and Mitohigh cells from wide-type regenerating planarians. The quantification is shown as mean ± SEM. n = 6 for 0 hpa and 4 for other time points. P-values were calculated by two-way ANOVA followed by Dunnett’s multiple comparisons.

To validate the enrichment of pluripotency in Mitolow cells in our study, we used the pan-neoblast marker piwi-1 to examine the piwi-1+ and piwi-1high cells within the five cell populations. A comparison among these populations showed that Mitolow cells contained the highest ratio of piwi-1+ cells (Fig. S1E, F) and enriched the highest proportion of piwi-1high neoblasts than the other cell populations as well as controls of Mitolow cells stained with piwi-1 sense probe or reagents without probes (Fig. 1D). We thereafter performed a transplantation assay to compare the pluripotency of Mitolow cells to Mitohigh cells. Following bulk-cell transplantation into neoblast-depleted planarian hosts, Mitolow cells proliferated more successfully than Mitohigh cells, forming colonies with an efficiency rate of 91% (n = 32). This rate surpassed that of the Mitohigh cells, which had a 4% efficiency rate (n = 25) (Fig. 1E). Additionally, only Mitolow cells rescued the viability of stem cell-depleted hosts that were irradiated with 60 Gy (Fig. 1F). Hence, these findings validated the pluripotency of Mitolow cells.

To explore the mitochondrial morphological dynamics in planarians, we examined whether mitochondria exhibit different morphologies in Mitolow, Mitohigh, CTGlow, and 2N cells. Leveraging the capacity of Hessian structured illumination microscopy (SIM), we directly observed the mitochondria in Mitolow, Mitohigh, CTGlow, and 2N cells. Our analysis revealed that Mitolow cells housed elongated, tubular mitochondria, whereas Mitohigh cells contained short, globular mitochondria (Fig. 1G). Quantitative analysis indicated a significant increase in mitochondrial size, total length, branch length, and ratio of length/branch of mitochondria from the Mitolow cells compared with the Mitohigh cells, consistent with the results of Hessian SIM (Fig. 1H–K). Mitochondrial morphology undergoes changes throughout the cell cycle. To further assess whether Mitolow and Mitohigh cells originate from similar cell cycles, we compared X1 cells from the S phase and G2/M phase with X2 and Xins cells. Our results revealed no obvious difference between cells in the S phase and G2/M phases, unlike those in X2 and Xins populations (Fig. S2A–E). This supports the aforementioned results of the cell cycle analysis, indicating that Mitolow and Mitohigh cells are in similar cell cycle stages. Additionally, we isolated Mitolow and Mitohigh cells from X1 to compare their mitochondrial morphology. We found similar mitochondrial morphology differences as we observed in the Mitolow and Mitohigh cells from SiRNeoblasts (Fig. S2F–J). By performing staining with MitoTracker Orange to measure the MMP, MitoTracker Green to measure the mitochondrial mass, and MitoSOX to measure reactive oxygen species (ROS) levels, we confirmed that Mitolow and Mitohigh cells exhibited lower and higher MMP (Fig. S2K), mitochondrial mass (Fig. S2L), and ROS levels (Fig. S2M), respectively.

We then examined the dynamics of the Mitolow and Mitohigh cell ratio during regeneration. Quantitative analysis of FACS data showed that the Mitolow cell ratio increased to its peak at 12 hours post-amputation (hpa), while the Mitohigh cell ratio decreased initially and then increased after 12 hpa, reaching its peak at 3 days post-amputation (dpa) (Fig. 1L, M). During planarian regeneration, neoblasts respond to injury by undergoing hyperproliferation, resulting in two peaks at 6 and 48 hpa, respectively38. The peak ratios of Mitolow and Mitohigh cells at 12 hpa and 3 dpa, respectively, may be the consequence of the substantial proliferation and differentiation of neoblasts. These observations prompted us to hypothesize that these variations in mitochondrial morphology might be related to changes in planarian cell fate status during body regeneration. Meanwhile, we acknowledge the limitation of this method for detecting the mitochondria in cells in vivo. Tomograph and other higher resolution in situ methods for mitochondria in planarians need to be developed to validate the dynamics in the future.

The mitochondrial fusion gene opa1 is required for planarian regeneration

To investigate the hypothesis suggested by our staining results, we next examined the relationship between mitochondria and regeneration, particularly the role of mitochondrial morphology changes during regeneration. The homologous gene candidates of six mitochondrial dynamics-related genes were identified in planarian S. med., including mitochondrial fusion genes Smedopa1, Smedmfn1, Smedmfn2, mitochondrial fission genes Smeddrp1 and Smedfis1, and mitochondrial transportation gene Smedmiro-1 (henceforth referred to as opa1, mfn1, mfn2, drp1, fis1, and miro-1) (Supplementary Data 1). Bacteria expressing dsRNA were used to feed planarians to knock down gene expression, and regeneration was tested after amputation. Through assessment for defects in regeneration but not homeostatic viability, we discovered that the mitochondrial fusion gene opa1 was essential for planarian regeneration postamputation, while knockdown of the mitochondrial fission gene drp1 resulted in only a minor regeneration delay for tail fragments (Fig. 2A, Supplementary Data 1). Quantitative real-time PCR (qPCR) showed a substantial reduction in transcripts after opa1 RNAi and drp1 RNAi (Fig. S3A, B). We then investigated organ integrity at 7 and 14 dpa via in situ hybridization and found that regeneration of the central nervous system (chat+), intestine (gata4/5/6+), and pharynx (laminin+) was blocked in planarians with opa1 knockdown compared with egfp RNAi controls (Fig. 2B, Fig. S3C).

Fig. 2: opa1 RNAi impairs planarian regeneration and induces alterations in mitochondrial dynamics.
figure 2

A Representative images show regeneration outcomes in live animals of egfp RNAi, opa1 RNAi, and drp1 RNAi animals at 7 dpa. Dpa, days post-amputation. n = 10 for each RNAi treatment. Scale bar, 500 μm. B In situ probes for chat and gata4/5/6 were used to show regeneration defects of the central nervous system and intestine in opa1 RNAi animals compared with egfp RNAi and drp1 RNAi animals at 7 dpa. n = 4 for each condition. Scale bar, 500 μm. C TEM images of planarian neoblasts (upper row) display the morphology of mitochondria in egfp RNAi, opa1 RNAi, and drp1 RNAi animals. Cells with neoblast characteristics were selected for observation. The experiment was repeated three times. The lower-row images correspond to the boxed areas in the upper-row images. Scale bars, 1 μm. D Representative Hessian-SIM images of mitochondria in Mitolow cells from egfp RNAi, opa1 RNAi, and drp1 RNAi planarians. Scale bar, 1 μm. EH Dot plots of mitochondrial size (E), total length (F), branch length (G), and length/branch ratio (H) of Mitolow cells from egfp RNAi (n = 24), opa1 RNAi (n = 38), and drp1 RNAi (n = 35) planarians. The quantification is shown as mean ± SEM. Unpaired two-tailed Student’s t test was used to determine the significance of differences between two conditions. I Pseudocolor dot plots show Mitolow and Mitohigh cells with their percentages from egfp RNAi and opa1 RNAi animals. Mitolow: p-value = 2.56 × 10-5 (72.9 ± 1.13% vs. 88.8 ± 0.80%); Mitohigh: p-value = 7.78 × 10-5 (20.6 ± 0.67% vs. 6.4 ± 0.74%) as opa1 RNAi versus egfp RNAi controls. The quantification is shown as mean ± SEM. Unpaired two-tailed Student’s t test was used to determine the significance of differences between two conditions.

To validate the effects of opa1 RNAi on mitochondrial morphology, we first showed that the ectopic expression of Myc-tagged planarian Opa1 could localize to the mitochondria in mammalian cells (Fig. S3D) due to the lack of an antibody against Opa1 and technical limitations of exogenous gene expression in planarian cells. Thereafter, transmission electronic microscopy (TEM) was used to compare the mitochondrial morphology in tissue sections from opa1 RNAi, drp1 RNAi, and egfp RNAi control animals. Previous studies have indicated that neoblasts could be identified with TEM through their lack of cytoplasmic differentiation and the presence of prominent decondensed chromatin30,39. Therefore, to compare among similar cell types, mitochondria were compared in neoblasts and were found that mitochondria within opa1 RNAi planarians were notably smaller and had lower cristae maturation than those within the egfp RNAi controls (Fig. 2C). In contrast, we observed increased mitochondrial interconnection and size in drp1 RNAi planarians compared with egfp RNAi controls (Fig. 2C). We also examined the consequence of mitochondrial morphology in mfn1;mfn2 RNAi, and fis1 RNAi animals. Despite significant reduction in their transcript levels (Fig. S3E, G), the mitochondrial morphology did not display any obvious changes in mfn1;mfn2 RNAi and fis1 RNAi compared to those in opa1 RNAi and drp1 RNAi animals (Fig. S3F, H). The lack of effect on mitochondrial morphology may explain the absence of defective regeneration in these RNAi conditions.

To confirm the changes in mitochondrial morphology elicited by opa1 RNAi and drp1 RNAi, we directly observed the mitochondria by Hessian SIM in isolated Mitolow cells from egfp RNAi, opa1 RNAi, and drp1 RNAi planarians. Our analysis revealed that egfp RNAi Mitolow cells housed elongated, tubular mitochondria, whereas opa1 RNAi Mitolow cells contained short, globular mitochondria (Fig. 2D) with a significant reduction in mitochondrial size (Fig. 2E), total length (Fig. 2F), branch length (Fig. 2G), and the ratio of length/branch of mitochondria (Fig. 2H). In contrast, drp1 RNAi Mitolow cells contained mitochondria with a significant increase in mitochondrial size, total length, and branch length compared with egfp RNAi Mitolow cells (Fig. 1D–H). These results confirmed the impact on the mitochondrial morphology caused by opa1 RNAi and drp1 RNAi. Concurrent with these alterations in mitochondrial dynamics, the application of opa1 RNAi significantly reduced the ratio of Mitolow cells compared to the egfp RNAi controls (Fig. 2I). These results suggest that defects in mitochondrial fusion by opa1 RNAi impact the pluripotency of Mitolow cells.

Reducing drp1 expression mitigates opa1 RNAi-induced regeneration impairment and mitochondrial morphological disruption

To examine whether mitochondrial dynamics plays a decisive role in cell fate and regeneration, we sought to ameliorate the deficiencies observed in animals subjected to opa1 RNAi. The drp1 RNAi resulted in an opposite morphology compared with that after opa1 RNAi (Fig. 2C–H). We, therefore, investigated whether drp1 RNAi ameliorated mitochondrial anomalies and regeneration defects in opa1 RNAi planarians. Subsequently, in addition to performing egfp;egfp RNAi, opa1;egfp RNAi, and drp1;egfp RNAi, we generated planarians with simultaneous opa1 and drp1 knockdown (referred to as opa1;drp1 RNAi animals). Unlike prior observations with opa1 RNAi animals, opa1;drp1 RNAi animals showed noticeable regeneration following amputation, mirroring the level of regeneration observed in the egfp;egfp RNAi and drp1;egfp RNAi animals (Fig. 3A). Probing for the intestinal and central nervous systems at 7 dpa via fluorescence in situ hybridization of mat and pc2 confirmed the rescue of the regenerative defects in opa1;drp1 RNAi animals (Fig. 3B). To validate the efficiency of gene suppression, we performed the qPCR analysis. The results showed that opa1;egfp RNAi and drp1;egfp RNAi significantly decreased the transcript levels of opa1 and drp1, respectively, compared to the egfp;egfp RNAi condition (Fig. S4A, B). Meanwhile, opa1;drp1 RNAi exhibited a comparable reduction in the transcript level of opa1 and drp1 to opa1;egfp RNAi and drp1;egfp RNAi conditions, respectively (Fig. S4A, B).

Fig. 3: Mitigation of opa1 RNAi-induced regeneration impairment and mitochondrial morphological disruption by reducing drp1 expression.
figure 3

A Images of live egfp;egfp RNAi, drp1;egfp RNAi, opa1;egfp RNAi, and opa1;drp1 RNAi animals at 7 dpa. n = 30 for each treatment. Scale bar, 500 μm. B In situ probes for mat and pc2 show the regeneration outcomes of intestinal, and neuronal cells, respectively, in egfp;egfp RNAi, opa1;egfp RNAi, and opa1;drp1 RNAi animals at 7 dpa. n = 3 for each treatment. Scale bar, 500 μm. C TEM images of planarian neoblasts (upper row) display the morphology of mitochondria in egfp;egfp RNAi, drp1;egfp RNAi, opa1;egfp RNAi, and opa1;drp1 RNAi animals. The lower-row images correspond to the boxed areas in the upper-row images. The experiment was repeated three times. Scale bars, 1 μm (upper row) and 200 nm (lower row). D Representative Hessian-SIM images of mitochondria in egfp;egfp RNAi, drp1;egfp RNAi, opa1;egfp RNAi, and opa1;drp1 RNAi Mitolow cells. Scale bars, 1 μm. EH Dot plots of mitochondrial size (E), total length (F), branch length (G), and length/branch ratio (H) egfp;egfp RNAi (n = 30), drp1;egfp RNAi (n = 28), opa1;egfp RNAi (n = 29), and opa1;drp1 RNAi (n = 35) Mitolow cells. The quantification is shown as mean ± SEM. Unpaired two-tailed Student’s t test was used to determine the significance of differences between two conditions.

We then investigated the phenotypic consequences of concurrent opa1 and drp1 knockdown on mitochondrial morphology via TEM. The results showed rescue of mitochondrial morphology in neoblasts of the opa1;drp1 RNAi animals similar to that in the egfp;egfp RNAi animals but in contrast to those in the opa1;egfp RNAi cohort (Fig. 3C). Furthermore, we directly observed mitochondria by Hessian SIM in isolated Mitolow cells from egfp;egfp RNAi, drp1;egfp RNAi, opa1;egfp RNAi, and opa1;drp1 RNAi planarians (Fig. 3D). Our analysis revealed that the mitochondrial morphology index values, which included the mitochondrial size (Fig. 3E), total length (Fig. 3F), branch length (Fig. 3G), and ratio of length/branch of mitochondria (Fig. 3H), were increased in the opa1;drp1 RNAi Mitolow cells to a similar level as in the egfp;egfp RNAi Mitolow cells, in contrast to those in the opa1;egfp RNAi Mitolow cells. A similar consequence was also observed in Mitohigh cells (Fig. S4C–G). The effect on the ratio of Mitolow and Mitohigh cells in different RNAi conditions was also assessed. Our analysis revealed that the ratio of Mitolow cells decreased concurrently with alterations in mitochondrial dynamics in both homeostasis (Fig. S4H, I) and at 1 dpa of opa1;egfp RNAi worms (Fig. S4J, K). This ratio was restored in opa1;drp1 RNAi worms. In contrast to the decrease in Mitolow cells, there was an increase in Mitohigh cells in opa1;egfp RNAi worms, which was also restored in opa1;drp1 RNAi Mitohigh cells.

Mitochondrial fusion–fission equilibrium is essential for neoblast proliferation and differentiation potency

To further evaluate the phenotypes of opa1 RNAi on regeneration, we examined its effects on cell proliferation and differentiation. In physiological homeostasis, no discernible decrease in neoblasts (Fig. 4A) or H3P+ proliferating neoblasts was detected upon opa1 RNAi (Fig. 4B). During regeneration in wild-type S. med., neoblasts are known to initiate two substantial proliferation bursts at intervals of 6 and 48 hpa38. We next examined the consequence of opa1 RNAi on neoblast proliferation during regeneration. In amputated opa1 RNAi worms, the ratios of H3P+ cells were significantly reduced at both 6 and 48 hpa compared to those in egfp RNAi controls (Fig. 4C–E), indicating a deficiency in the capacity of neoblasts to respond to tissue loss. In line with this restoration of mitochondrial morphology, cell proliferation rates during regeneration were also increased in opa1;drp1 RNAi animals compared with opa1;egfp RNAi animals according to the density of H3P+ cells at 6 and 48 hpa, respectively (Fig. 4D, E).

Fig. 4: Neoblast proliferation and differentiation potency is compromised by opa1 RNAi.
figure 4

A Representative images of FISH staining for piwi-1 in egfp RNAi and opa1 RNAi planarians in homeostasis at 7 days post-feeding (dpf). Scale bar, 200 μm. B Images of egfp RNAi and opa1 RNAi planarians at 7 dpf stained with an anti-H3P antibody. Scale bar, 200 μm. C Changes in H3P+ cell density during regeneration of different RNAi planarian tails. n = 4 tails biological replicates from each RNAi condition at every time point. D, E Changes in H3P+ cell density at 6 hpa (D) and 48 hpa (E) for the indicated RNAi planarians. Besides drp1 RNAi animal tail samples at 48 hpa (n = 4), n = 5 tail samples were used from each of the other RNAi conditions at 6 hpa and 48 hpa. Each dot indicates an independent individual animal. F Representative images of piwi-1+ signals in different RNAi planarians at 7 and 14 dpi following 12.5 Gy of irradiation. Scale bars, 500 μm. G Quantification of piwi-1+ cell density at 14 dpi following 12.5 Gy of irradiation of the indicated RNAi animals. Besides opa1 RNAi animal samples (n = 9), n = 10 independent individual samples from each of the other RNAi conditions. HJ Quantification of co-stained cells positive for PIWI-1 and ston2 (H), hnf4 (I), or prog-1 (J) at 2 dpa in different RNAi planarians. Besides opa1 RNAi animal samples (n = 4), n = 3 samples from each of the other RNAi conditions (H), and n = 4 biological replicates from each RNAi condition (IJ), each dot indicates an independent data point. K Bar plot shows the ratio of X1 and X2 cells in different RNAi planarians. n = 4 biological replicates for each RNAi condition. L Bar plot shows the ratio of S and G2/M stages in X1 cells from different RNAi planarians. n = 4 biological replicates for each RNAi condition. The quantification in the analysis of all data in Fig. 4 (CE) and (GL) is shown as mean ± SEM, a two-tailed Student’s t-test was used to calculate the p-value versus egfp RNAi controls.

Consistently, we noticed significant suppression of neoblast repopulation following sublethal irradiation in organisms subjected to opa1 RNAi compared to egfp RNAi controls (Fig. 4F, G; Fig. S5A, B). Neoblast repopulation was increased to normal in opa1;drp1 RNAi animals compared with egfp;egfp RNAi and drp1;egfp RNAi controls (Fig. 4F, G). To determine whether the defects in neoblast proliferation are cell autonomous, we transplanted neoblasts isolated from opa1 RNAi planarians into neoblast-depleted wild-type hosts and vice versa and performed the same procedure. As predicted, the neoblasts of opa1 RNAi planarians failed to exhibit engraftment and clonogenic capacity in contrast to neoblasts of egfp RNAi control animals, as indicated by the significant difference in piwi-1+ cell number at 7 days post-transplantation (7 dpt) (Fig. S5C, D). In contrast, transplanted healthy, wild-type neoblasts engrafted and proliferated in irradiated opa1 RNAi worms at a level higher than that in irradiated egfp RNAi worms (Fig. S5E). These results suggest that Opa1-regulated mitochondrial fusion is cell-autonomously essential for maintaining neoblast proliferation in response to regeneration.

Interestingly, when we examined whether drp1 RNAi-induced mitochondrial interconnection affects neoblast proliferation during regeneration, we found significant reductions in H3P+ cells at both 6 and 48 hpa (Fig. 4D, E), but only slightly reduced repopulation after sublethal irradiation (Fig. 4F, G; Fig. S5A, B) and no difference after transplantation of drp1 RNAi neoblasts compared to egfp RNAi controls (Fig. S5C–E). Considering the successful regeneration of drp1 RNAi planarians, these results suggest that the reduced proliferation by drp1 RNAi is a non-cell-autonomous effect on neoblasts and that reduced proliferation alone does not explain the regeneration defects in opa1 RNAi planarians.

To further elucidate the function of mitochondria in cell differentiation during regeneration, we harvested regenerating trunks at 2 dpa and proceeded with whole-mount labeling of the neural progenitor marker ston2, intestinal progenitor marker hnf4, and early epidermal progenitor marker prog1. Immunofluorescence staining revealed significantly lower levels of ston2, hnf4, and prog-1 co-expression with PIWI-1 in the regenerating blastema region upon opa1 RNAi compared with egfp RNAi controls (Fig. 4H–J, Fig. S5F–H). The differentiation rates during regeneration were also increased in opa1;drp1 RNAi animals compared with opa1;egfp RNAi animals according to the density of progenitors of ston2+PIWI-1+ (Fig. 4H), hnf4+PIWI-1+ (Fig. 4I), and prog-1+PIWI-1+ (Fig. 4J) cells in the blastema at 2 dpa. Staining the patterning markers at 3 dpa for anterior sfrp-1 and notum and posterior wnt-1, wnt11-2, and wntP-2 revealed that the reestablishment of tissue patterning was disturbed in opa1 RNAi planarians compared to egfp RNAi controls (Fig. S5I). To assess whether the defects in cell differentiation could be attributed to an influence on the cell cycle, we performed cell cycle analysis on the proportion of X1 cells and the S/G2/M stages. The statistical results showed no significant reduction in neoblasts in opa1;egfp RNAi planarians on the proportion of X1 and X2 cells (Fig. 4K, Fig. S5J). Further analysis of the cell cycle of X1 cells revealed a decrease in the ratio of cells in the G2/M phase and an increase in the S phase (Fig. 4L, Fig. S5K), which potentially influenced the cell fate in opa1;egfp RNAi planarians. These above findings indicate that maintaining the equilibrium of mitochondrial dynamics is fundamental for determining stem cell fate when these cells must replace missing tissues and meet the demands of regeneration.

Mitolow cells are enriched at earlier stages of fate choices than Mitohigh Cells

Before delving into the role of opa1 in neoblast pluripotency, we first conducted experiments to elucidate the relationship between neoblast pluripotency and mitochondrial dynamics. After isolating cells from wild-type intact planarians using flow cytometry, we employed single-cell RNA sequencing (scRNA-seq) and subjected the cells to UMAP analysis, which resulted in 23 distinct clusters (Fig. 5A). Each cluster was assigned to a specific cell type based on the expression of neoblast subtype markers40,41,42,43,44,45 (Fig. 5A, Fig. S6A, Supplementary Data 2). Subsequently, we determined the lineage identity of each cell and compared the compositions of Mitolow and Mitohigh cells. Within the identified clusters, two clusters showed high expression of tgs-1, which we referred to as tgs-1+ neoblast c2 and c3 for further analysis. Analysis of individual scRNA-seq datasets revealed distinct cluster compositions between Mitolow and Mitohigh cells (Fig. 5B). While both Mitolow and Mitohigh cells comprised epidermal, muscular, cathepsin+/phagocytotic, parenchymal, protonephridia, and neural neoblasts, they differed in the composition of these clusters. Mitolow cells predominantly contained tgs-1+ neoblast c2 and c3, as well as intestinal neoblasts. In addition, Mitolow cells expressed higher levels of neoblast-enriched genes and lower levels of XBP1 and P4HB compared with Mitohigh cells (Fig. S6B), which matches the gene expression signatures of pluripotent neoblasts in the single-step fate model46. These results provide evidence to support the higher pluripotency status of Mitolow cells. To avoid the batch effect, we performed additional scRNA-seq analysis on Mitolow and Mitohigh cells isolated from tails and intact animals. The results showed a similar distinguishment in the enrichment of different cell clusters (Fig. S7A).

Fig. 5: Mitolow cells are enriched with neoblasts at earlier stages of fate choices than Mitohigh cells.
figure 5

A UMAP plot of integrated single-cell data from Mitolow and Mitohigh cells. The cell types of each cluster in the same lineage are grouped with dashed shapes. B UMAP plots show clusters enriched in Mitolow and Mitohigh cells, respectively. The clusters grouped into the same cell lineages are shown in the same color. C UMAP plots show the subclustered epidermal neoblasts in Mitolow and Mitohigh cells, respectively. D Dot plot shows the expression of epidermal subtype markers across epidermal Mitolow and Mitohigh clusters. E Expression of subtype-specific markers in epidermal Mitolow and Mitohigh clusters. F UMAP plots show the subclustered muscular neoblasts in Mitolow and Mitohigh cells, respectively. G Dot plot shows the expression of muscular subtype markers across muscular Mitolow and Mitohigh clusters. H UMAP plots show the subclustered cathepsin+ neoblasts in Mitolow and Mitohigh cells, respectively. I Dot plot shows the expression of cathepsin+ subtype markers across cathepsin+ Mitolow and Mitohigh clusters. J UMAP plot shows the subclustered tgs-1+ neoblast cluster 2 (c2). K Expression of subtype-specific markers in the subclustered tgs-1+ neoblast c2. L UMAP plot shows the subclustered tgs-1+ neoblast c3. M Expression of subtype-specific markers in the subclustered tgs-1+ neoblast c3.

We next sought to elucidate the transcriptomic signatures associated with the pluripotency of Mitolow cells. We utilized the recently reported fate-specific transcription factors (FSTFs) to validate the subclusters within the epidermal, muscular, cathepsin+, and neural lineages33. By integrating and subclustering the epidermal cells from the Mitolow and Mitohigh cells, we identified 11 subclusters of epidermal neoblasts (Fig. 5C). The Mitolow epidermal neoblasts exhibited higher expression of canonical FSTFs such as zfp-1, co-expressing with soxP-3, and p53, while persisting the expression of soxP-3, and p53 in Mitohigh epidermal neoblasts (Fig. 5D, E). Conversely, the Mitohigh epidermal neoblasts predominantly expressed G0 progenitor FSTFs such as prog-2 and agat-3 (Fig. 5D, E). Notably, the dorsal and ventral epidermis, marked by the expression of PRDM-1 and kal1 respectively, were enriched in Mitohigh cells (Fig. 5D, E). These results suggest that the Mitolow epidermal neoblasts are more enriched at earlier stages of fate choices compared to Mitohigh cells. Subclustering muscular neoblasts revealed a total of 18 subclusters (Fig. 5F). Clusters 0, 1, 2, 3, 5, and 15 were enriched in Mitolow cells and displayed highly expressed HIST1H2AA and H2AFJ genes, which are core components of the nucleosome. In contrast, most subclusters in Mitohigh cells exhibited enriched expression of cell-type-specific FSTFs, including longitudinal muscle marker myoD, circular muscle marker nkx1-1, and intestinal muscle marker gata4/5/6-3 and PTPRD. Importantly, these TFs were already expressed during the Mitolow stage, suggesting the onset of muscular neoblast diversity within Mitolow cells (Fig. 5G, Fig. S7B). Our subclustering of cathepsin+ neoblasts revealed 5 clusters (Fig. 5H). Cluster 0, enriched in Mitolow cells, exhibited high expression of zfp-1 and TYMS genes. Conversely, the cathepsin+ cell marker CTSL was enriched in Mitohigh neoblasts, suggesting the enrichment of earlier phagocytic cell fates in Mitolow cells (Fig. 5I, Fig. S7C). When examining the subclusters of tgs-1+ neoblast c2 and tgs-1+ neoblast c3 within Mitolow cells, we identified 8 subclusters (Fig. 5J, L). Compared to tgs-1+ neoblast c3, tgs-1+ neoblast c2 had a higher enrichment of tgs-1+, gcm-1+, and ski-3+ cells. Of particular note, the neural FSTFs such as pou4l-1 and neuroD-1 were highly expressed in tgs-1+ neoblast c2 (Fig. 5K, M). Earlier stage FSTFs, such as the epidermal marker zfp-1, muscular marker myoD, and intestinal marker gata4/5/6-2, were enriched in the tgs-1+ neoblast c3 subcluster (Fig. 5K, M). Previous research proposed that tgs-1+ neoblasts might be neural-specialized neoblasts46. Our dataset, however, revealed previously undescribed cell clusters that co-express tgs-1 with multiple FSTFs. The lineage hierarchy of these cells warrants further detailed analysis in future studies. Altogether, our cell atlas of Mitolow and Mitohigh cells indicates that Mitolow cells show higher enrichment at earlier stages of fate choices compared to Mitohigh cells, as evidenced by the differential expression of FSTFs across various lineages.

Changes in mitochondrial dynamics regulate the gene expression in multiple lineages of neoblasts

To investigate gene expression changes in RNAi conditions, we performed bulk RNA-seq and scRNA-seq analyses (Fig. 6A). Bulk RNA-seq was conducted on egfp RNAi and opa1 RNAi tripartite fragments at 1 dpa as well as on tail fragments at various time points post-amputation (Fig. 6A). The results of bulk RNA-seq analysis were consistent with qPCR, indicating that the double knockdown of opa1 and drp1 effectively reduced the transcript levels of both genes to a similar extent as their respective single knockdown conditions (Fig. S8A, B). To assess the wound response and neoblast response in the early stages of regeneration, we analyzed the expression of wound-induced genes and neoblast-enriched genes. The bulk RNA-seq analysis showed that opa1 RNAi animals, which can normally activate the expression of early wound-induced genes, exhibited higher expression levels of pim-3, runt-1, and wnt-1 at 1 dpa. (Fig. 6B, Supplementary Data 3). The heatmap also showed a decrease in the expression levels of neoblast-enriched genes, particularly piwi-1 and histone2B, in opa1 RNAi animals compared to egfp RNAi animals at 1 dpa. (Fig. 6B, Supplementary Data 3). Additionally, piwi-1 and histone2B exhibited reduced transcript levels at 2 dpa, along with other neoblast-enriched transcripts including mcm4, vasa-1, and cyclinB-1 (Fig. S8C, Supplementary Data 3). These results are consistent with the findings from H3P antibody staining experiments conducted at both 6 and 48 hpa in opa1 RNAi animals compared to egfp RNAi controls (Fig. 4C–E). It further supports the result that there is a significant deficit in neoblasts’ capacity to respond to injury. Furthermore, the bulk RNA-seq data showed a decrease in the expression of anterior PCGs, including ndl-2, wnt2, ndl-5, fz5/8-4, sfrp-1, ndk, and fz4-4 in opa1 RNAi (Fig. S8C, Supplementary Data 3), consistent with the FISH results that opa1 RNAi animals were unable to establish the positional information during regeneration (Fig. S5I).

Fig. 6: Modulations in mitochondrial dynamics influence gene expression across diverse neoblast lineages.
figure 6

A Illustration of gene expression analysis for bulk RNA-seq of tripartite fragments at 1 dpa from egfp RNAi, opa1 RNAi, and drp1 RNAi, tail fragments at various time points post-amputation from egfp RNAi and opa1 RNAi, and Mitolow scRNA-seq at 1 dpa. B Heatmap shows the expression of signature genes after opa1 knockdown. The red asterisk (*) means DEGs compared opa1 RNAi with egfp RNAi control group. The up- (rose boxes) and down- (black boxes) regulated genes were identified with the threshold of Fold change > 1.2 and adjusted p-value < 0.05. C UMAP plot of integrated single-cell data from Mitolow cells of the indicated RNAi animals at 1 dpa. The clusters categorized into the same cell types are shown in the same color. D Differential gene expression analysis comparing opa1;egfp RNAi and egfp;egfp RNAi Mitolow cells shows up- and down-regulated genes across all cell clusters. Top: Bar graph showing the number of DEGs in each cluster. Middle: For each cluster, the change in expression of significantly changed genes (p adj <0.05, |log2FC | > 0.585) is shown. Bottom: The cell fate to which different clusters belong. The colors indicate the clusters. LogFC, log fold change. E Venn diagram presents the number of genes significantly downregulated in opa1;egfp RNAi compared with egfp;egfp RNAi, but significantly upregulated in opa1;drp1 RNAi compared with opa1;egfp RNAi Mitolow cells. F KEGG pathway enrichment and (G) Protein-protein interaction networks (PPI) display genes selected from the Venn diagram in (E). The Benjamini-Hochberg method was used to adjust the p values. The adjusted p-value < 0.05 was set to determine the significantly changed expression level. H Venn diagram presents the number of genes significantly upregulated in opa1;egfp RNAi compared with egfp;egfp RNAi, but significantly downregulated in opa1;drp1 RNAi compared with opa1;egfp RNAi Mitolow cells. I KEGG pathway enrichment and (J) Protein-protein interaction networks (PPI) display genes selected from the Venn diagram in (H). The Benjamini-Hochberg method was used to adjust the p values. The adjusted p-value < 0.05 was set to determine the significantly changed expression level.

To investigate the impact of opa1 knockdown on neoblast pluripotency and the proliferation reduction observed during regeneration, we compared the scRNA-seq data of Mitolow cells at 1 dpa from animals subjected to egfp;egfp RNAi, drp1;egfp RNAi, opa1;egfp RNAi, and opa1;drp1 RNAi. In total, 39,596 cells were collected from all four groups at 1 dpa, and egfp;egfp RNAi were performed twice as a technical replicate. Cells were categorized into 17 clusters by utilizing UMAP analysis on the expression of enriched genes (Fig. 6C). The same gene signatures that were employed to annotate the Mitolow cells from wild-type intact planarians were used to confirm the identity of each cluster (Fig. S9A, Supplementary Data 2). Consistency in cell clusters and cell ratios between the two egfp RNAi replicates indicated the robustness of the experimental results (Fig. S9B, C). Compared to Mitolow cells from wild-type intact planarians, the integrated cell atlas from the above RNAi animals at 1 dpa also comprised epidermal, protonephridia, intestinal, cathepsin+, muscular, parenchymal, and tgs-1+ neoblasts. Notably, the extended source of cells enabled the distinction of foxA+ and NKX2-3+ parenchymal neoblasts (Fig. 6C). We also identified a previously unreported cell cluster characterized by the expression of runt-1. Since runt-1 is a wound-induced gene in neoblasts47, this cell population might represent a specific neoblast population induced in regeneration. Furthermore, we also identified an unknown neoblast cluster exhibiting high expression of histone-encoding genes HIST1H2AA and H2AFJ. These newly identified cell clusters present intriguing targets for future studies into their functions in planarians.

To understand the impact of RNAi treatment, we compared the proportion of each cell cluster across different RNAi groups. No obvious difference was observed (Fig. S9B, C). We next analyzed the differential expressed genes (DEGs) in each cell cluster by comparing them with egfp;egfp RNAi controls (Supplementary Data 4). We identified 152 upregulated genes and 760 downregulated genes when comparing opa1;egfp RNAi with egfp;egfp RNAi (Fig. 6D). A noticeable fewer DEGs were observed in drp1;egfp RNAi (Fig. S9D), supporting the findings of minor influence on neoblast pluripotency by drp1 knockdown (Fig. 3D–G). To further explore the rescue of gene expression, we compared the gene expression in opa1;drp1 RNAi with opa1;egfp RNAi (Supplementary Data 4). 61 genes that were downregulated in opa1;egfp RNAi were significantly increased by opa1;drp1 RNAi (Fig. 6E), including genes associated with oxidative phosphorylation (OXPHOS), mitochondrial biogenesis, and PPAR signaling pathways. Importantly, the expression of 60 out of these 61 genes was rescued in opa1;drp1 RNAi by comparing with egfp;egfp RNAi treatment (Supplementary Data 4). The recovery of genes enriched in the ribosome and PPAR signaling pathway suggests that the cellular translation activity was restored (Fig. 6F, G). Intriguingly, egr-4, which is essential for planarian head regeneration48, was downregulated in opa1;egfp RNAi but rescued in opa1;drp1 RNAi. Moreover, fos-1 as one of the first wave genes in response to wound47, was downregulated in opa1;egfp RNAi and restored in opa1;drp1 RNAi planarians. This suggests that the mitochondrial dynamics, affected by opa1 and drp1, are also involved in the regulation of early response signals that are critical for planarian regeneration. 59 genes shown as upregulated in opa1;egfp RNAi were significantly decreased by opa1;drp1 RNAi (Fig. 6H). These genes, including RARS1, NARS1, DARS1, MARS1, and EPRS1 in the aminoacyl-tRNA biosynthesis and CARS in the transfer RNA biogenesis pathways, were enriched in these DEGs (Fig. 6I, J; Supplementary Data 4). A group of the upregulated DEGs (19/59) were rescued in cluster 11 muscular neoblasts, like RPS29, a 40S ribosome protein. This suggests that the opa1;egfp RNAi causes abnormal translation, leading to upregulated transcription of genes associated with the tRNA biosynthesis, which is restored in opa1;drp1 RNAi planarians.

To assess the impact of opa1 RNAi across all cell types besides stem cells, we compared the single-cell atlas of egfp;egfp RNAi, drp1;egfp RNAi, opa1;egfp RNAi, and opa1;drp1 RNAi animals. Given the potential of dissociation manipulation to result in the loss of fragile cells and induce stress responses, we employed single-nucleus RNA sequencing (snRNA-seq) to characterize the cellular components in RNAi planarians. This analysis yielded a total of 42,965 nuclei, which were subsequently categorized into 26 clusters through UMAP based on the expression of enriched genes in each cluster. The average number of genes per cell nucleus was 1,031, consistent with the typical range for snRNA-seq. Using marker genes, we identified 9 major cell types in the cell atlas (Fig. S10A, B, Supplementary Data 2). These clusters were then validated using classic planarian cell type markers (Fig. S10A), including neoblast, epidermal, muscular, intestinal, cathepsin+/phagocytic, protonephridia, pharyngeal, neural, and parenchymal cells. The diversity of cell types and gene heterogeneity were consistent with previous reports42,43.

We also compared the proportion of each cell cluster among different RNAi groups. No obvious difference was observed (Fig. S10C). We then examined the DEGs in each cell cluster by comparing them with egfp;egfp RNAi controls (Fig. S11A–C, Supplementary Data 5). We identified 10 downregulated genes in opa1;egfp RNAi that were rescued by opa1;drp1 RNAi (Fig. S11D, E). These included MT-ND1, MT-ND3, ND4, ATP6, MT-COX1, MT-COX2, and MT-COX3, which are related to OXPHOS and encoded by the mitochondrial genome (Supplementary Data 5). These genes were downregulated in epidermal, neoblast, and intestinal clusters, but their expression was restored in neoblasts following opa1;drp1 RNAi (Fig. S11A, D, E; Supplementary Data 5). Similarly, we observed these genes downregulated in 1 dpa mitolow scRNA-seq analysis after opa1 knockdown (Fig. 6E, Supplementary Data 4). Furthermore, we identified 12 upregulated genes that were rescued by opa1;drp1 RNAi (Fig. S11F, Supplementary Data 5). Genes encoding IARS1, QARS1, AARS1, and RARS1, which are involved in aminoacyl-tRNA biosynthesis, were upregulated in neoblasts, epidermal cells, and intestinal cells after opa1 knockdown. Their expression was restored in epidermal cells following opa1;drp1 RNAi. This indicates that opa1 knockdown stimulates these genes specifically in epidermal cells. Despite these findings, it is important to note that the number of genes differentially expressed following opa1 knockdown in other cell types was smaller than in neoblasts. This suggests that opa1-regulated mitochondrial dynamics primarily influence genes and functions involved in neoblasts during the early stages of regeneration, although we cannot exclude the plausibility of other cell types.

To complement the single-cell transcriptional analysis and identify gene expression changes specific to opa1 RNAi at 1 dpa, we further utilized bulk RNA-seq data from intact and 1 dpa time points (Fig. 6A). By comparing with egfp RNAi and drp1 RNAi, we identified 111 genes that were specifically downregulated, and 26 genes upregulated in opa1 RNAi planarians at 1 dpa (Fig. S12A–F, Supplementary Data 6). Pathway enrichment analysis found that UQCR10, UQCRH, and UQCRQ, which are associated with OXPHOS, were detected in bulk RNA-seq data (Fig. S12G, H). Histone proteins such as H4C1, H3-5, H2AZ1, and H3C1 showed lower expression in opa1 RNAi worms (Supplementary Data 6). In addition, genes involved in RNA splicing, including H4C1, YBX3, SNRPF, SNRPD1, SNRPD2, and SF3B5, were downregulated (Fig. S12H). The translation initiation factor EIF1, along with downregulated ribosomal protein genes RPLP2 and RPS25, suggested that opa1 knockdown decreases transcription and translation during regeneration (Fig. S12H). Furthermore, RARS1 associated with arginyl-tRNA was upregulated at 1 dpa (Fig. S12I), suggesting a feedback regulation to compensate for reduced translation following opa1 RNAi. These DEGs were enriched in similar signaling pathways compared to those enriched in the scRNA-seq analysis (Fig. 6F, G, I, J; Fig. S11E, F), offering potential insights into the regulatory pathways affected by mitochondrial fusion defects.

Mitochondrial Fusion–Fission Equilibrium is Critical for Mitonuclear Balance

To facilitate the study of DEGs related to mitochondrial dynamic changes, we comprehensively categorized the annotated genes encoded by the nuclear genome and mitochondrial DNA (mtDNA). Previous studies have estimated that the mitochondrial proteome in yeast and humans contains at least 1000  1500 different proteins49,50. In a study of the planarian mtDNA, it was found that the sexual and asexual biotypes of Schmidtea mediterranea have mtDNA with lengths of 27,133 bp and 17,176 bp, respectively. However, both biotypes have the same number of genes encoding for 12 proteins (ND1, ND2, ND3, ND4, ND4L, ND5, ND6, COX1, COX2, COX3, CYTB, and ATP6), 2 rRNA, and 22 tRNAs51. Although the planarian mitochondrial proteome has not been reported, gene homologs have been identified by searching the nuclear genome sequence for known mitochondrial-related genes (Supplementary Data 7)52. The identified proteins were then grouped according to their major functions, including mitochondrial detoxification, mitochondrial dynamics and surveillance, mitochondrial ribosome assembly, translation, OXPHOS, and mitochondrial metabolism. Consequently, we integrated and classified the genes that were transcriptionally up-regulated and down-regulated in association with opa1 RNAi, compared with egfp RNAi controls (Supplementary Data 7, Fig. 7A).

Fig. 7: Regulation of the mitonuclear balance by the equilibrium of mitochondrial dynamics is essential for regeneration.
figure 7

A Summary of the DEGs in multi-transcriptomic analysis. Genes are summarized with subcellular localization and annotated functions. Downregulated and upregulated genes are in blue and red, respectively. B Left: representative WB images of egfp;egfp RNAi, opa1;egfp RNAi, drp1;egfp RNAi, and opa1;drp1 RNAi planarians treated with puromycin (puro) or puromycin and cycloheximide (CHX). Ponceau staining was used to confirm equal protein loading. n  =  3 biological replicates for each condition; each biological replicate was the mixture of proteins of six RNAi worms. Right: Quantitative representation of puromycin incorporation observed in left WB images. The data are shown as the mean ± SEM. P-value was calculated by unpaired two-tailed Student’s t-test versus egfp;egfp RNAi controls. C Quantification of the relative expression level of UPR signals in opa1 RNAi animals compared with egfp RNAi animals. n = 3 biological replicates for each RNAi condition. The quantification is shown as mean ± SEM. P-value was calculated by unpaired two-tailed Student’s t-test versus egfp RNAi controls. D Representative live images show the block of regeneration after mrpl2 RNAi and uqcr10 RNAi. Scale bar, 500 μm. E Quantification of H3P+ cell density at indicated regeneration time points of egfp RNAi, mrpl2 RNAi, and uqcr10 RNAi animals. Besides uqcr10 RNAi animal samples at 48 hpa and 72 hpa (n = 4 biological replicates for each condition), n = 3 samples from each of the other RNAi conditions at each time point. The data are shown as the mean ± SEM. Unpaired two-tailed Student’s t-test was used to determine the significance of differences at each time point versus egfp RNAi controls. F Representative FISH images of anterior and posterior patterning markers at 72 hpa. n = 3 worms of biological replicates for each RNAi treatment. Scale bar, 50 μm.

Through above transcriptome analysis, we found that the defects in opa1-regulated mitochondrial dynamics caused the downregulation of nuclear genome encoded H4C1, YBX3, SF3B5, SNRPD1, SNRPD2, and SNRPF, which are enriched in the U2-type prespliceosome assembly pathway, upregulation of nuclear genome encoded RARS1, EPRS1, IARS1, AARS1, and QARS1, which are enriched in the cytosolic tRNA aminoacylation signaling pathway, but also caused the downregulation of mitochondrial DNA encoded MT-ND1, MT-CO1, and ATP6 genes, which are involved in the electron transport chain, and MRPL2, MRPL17, and MRPL24 genes, which are involved in mitochondrial protein translation (Fig. 7A). Mitonuclear balance refers to the delicate equilibrium between the expression and function of genes encoded in the nuclear DNA and those encoded in the mitochondrial DNA (mtDNA)20,53,54. The above analyses suggest that there is an imbalance in the transcription of genes for mitochondrial proteins and cytosolic proteins (Fig. 7A). To assess the translation level, a puromycin incorporation experiment was performed. The translation inhibitor CHX was used as a negative control, and no significant decrease in puromycin incorporation was observed when comparing opa1;egfp RNAi with egfp;egfp RNAi control animals at 1 dpa (Fig. 7B). This suggests that there is no significant decrease in protein translation levels at the bulk-tissue level between opa1;egfp RNAi and egfp;egfp RNAi controls (Fig. 7B). However, considering the gene expression difference in neoblasts, we could not rule out the possibility that a reduction in protein translation occurred in neoblasts and other specific cell types.

Proper mitochondrial fusion and fission play pivotal roles in maintaining mitonuclear balance, facilitating coordinated gene expression, efficient energy production, and effective stress responses9,55. Previous studies have shown that loss-of-function of opa1 could result in the downregulation of mitochondrial genes encoding for components of the OXPHOS, mitochondrial ribosomal proteins, and other mitochondrial biogenesis-related genes, and upregulation of genes of heat shock proteins (HSPs), reactive oxygen species (ROS) scavengers, pro-apoptotic factors, and genes involved in UPR56,57,58,59,60,61. Knockout of either opa1 or drp1 would result in an enhanced UPR55,62. To validate this response, we attempted to identify homologous genes involved in mito-nuclear signaling pathways, such as atfs-1, atf4/5, and Nrf2, that are found in the roundworm C. elegans and mice. However, no obvious homologs were identified, likely due to the divergence of the DNA sequences with evolutionary distance. (Supplementary Data 7)63. Instead, 29 other homologous genes associated with UPRmt and UPRER were identified in planarians64,65. These genes included HSP10, HSP60, eIF2α, sirt1, sirt3, catalase, PERK, and among others. This divergence can obscure the identification of direct orthologs. Despite this challenge, we were able to quantify the transcription levels of genes associated with the UPRER under opa1 RNAi condition (Fig. 7C). This approach not only validates our response but also underscores the complexity and conservation of mito-nuclear signaling pathways across different organisms. The lack of many core conserved components of UPRmt in planarians suggests that UPRmt regulation in these organisms may involve other, yet undiscovered, transcription factors.

Given the altered expression of nuclear genome-encoded mitochondrial subunits, we hypothesize that the changes in mitochondrial dynamics may decrease the function of the electron transport chain and mitochondrial ribosomal subunits, thereby exerting functional effects. To further validate our hypothesis, we then performed an RNAi screen targeting the DEGs identified from opa1 RNAi animals. We focused on the nuclear genome-encoded genes with homologs in mammals and ultimately identified 62 genes that met these criteria (Supplementary Data 8). RNAi knockdown of two genes, mrpl2 and uqcr10, which encode proteins located in mitochondria, resulted in defective regeneration without noticeable homeostasis defects (Fig. 7D, Supplementary Data 8). Furthermore, neoblast hyperproliferation during regeneration was reduced in mrpl2 RNAi and uqcr10 RNAi animals compared with egfp RNAi controls (Fig. 7E). Additionally, the expression patterns of sfrp-1, notum, and wnt-1 were disorganized (Fig. 7F), similar to the observations in opa1 RNAi animals. These results indicate that the genes regulated by mitochondrial dynamics indeed contribute to the function of tissue regeneration. Moreover, upregulated genes such as IARS1, QARS1, and RARS1 are also essential for planarian normal physiology (Supplementary Data 8) and likely act as compensatory mechanisms in response to changes in mitochondrial dynamics by increasing their expression levels. Taken together, our findings suggest a mechanism in which the mitochondrial morphology dynamics are crucial for coordinating the mitonuclear balance, likely through the regulation of metabolism and translation in mitochondria.

Mitochondrial fusion–fission equilibrium is critical for metabolism in regeneration

Consistent with previous studies60,66,67,68,69,70,71, our results showed the downregulation of mitochondrial genes encoding for components of the OXPHOS in opa1;egfp RNAi animals. To examine the influence on mitochondria, Seahorse assays were employed to measure the oxygen consumption rate (OCR) and extracellular acidification rate (ECAR). OCR measurements suggested a delayed response but comparable OXPHOS in opa1;egfp RNAi animals, which was restored in opa1;drp1 RNAi animals (Fig. S12J). ECAR measurements suggested an increase in glycolysis in opa1;egfp RNAi and drp1;egfp RNAi, but this was recovered in opa1;drp1 RNAi animals (Fig. S12K). ATP levels were measured and found to be significantly reduced in opa1;egfp RNAi animals at 1 dpa (Fig. 8A), and in all drp1;egfp RNAi, opa1;egfp RNAi, and opa1;drp1 RNAi animals compared to egfp RNAi controls at 3 dpa (Fig. 8B). These observations suggest that compromised mitochondrial dynamics reduced the function of OXPHOS in planarians.

Fig. 8: Dysfunction of metabolism after influencing the mitochondrial dynamics in planarians.
figure 8

A, B The bar plot shows the relative level of ATP in different RNAi planarians at 1 dpa (A) and 3 dpa (B). n = 4 independent biological replicates for each RNAi condition. The data are shown as the mean ± SEM. P-value was calculated by unpaired two-tailed Student’s t-test versus egfp;egfp RNAi controls. C Oil Red O staining for lipid droplets on the cryosections of the indicated RNAi planarians. Scale bar, 50 μm. D Bar plot shows the density of the lipid droplets in size and the ratio among three ranges of size ( < 6, 6–20, and > 20 μm2) in egfp;egfp RNAi (n = 4580), opa1;egfp RNAi (n = 5134), drp1;egfp RNAi (n = 6706), and opa1;drp1 RNAi (n = 4667). the n value represents the number of lipid droplets. Each dot indicates a lipid droplet. Data are shown as mean ± SEM P values were calculated by one-way ANOVA followed by Tukey’s multiple comparisons. E Body length of egfp;egfp RNAi (n = 17), opa1;egfp RNAi (n = 15), drp1;egfp RNAi (n = 15), and opa1;drp1 RNAi planarians (n = 17) at 7 dpa. The quantification is shown as mean ± SEM. P-value was calculated by unpaired two-tailed Student’s t-test. F Heatmap shows the changes of metabolites in indicated RNAi animals at 1 dpa. Four replicates for each RNAi were grouped in the indicated colors. The metabolites with significant changes in abundance (z score) in opa1;egfp RNAi compared to egfp;egfp RNAi were shown. The Benjamini-Hochberg method was used to adjust the p values. The adjusted p-value < 0.05 and fold change > 1.2 was set to determine the significantly changed expression level. G Cartoon illustration for different morphologies of mitochondria and pluripotency in Mitolow cells and Mitohigh cells. H Illustration for the regeneration defects by opa1 RNAi and rescue by opa1;drp1 RNAi. I The balance of mitochondrial fusion and fission is regulated by opa1 and drp1. J opa1 RNAi causes the malfunction of mitochondria, elevated UPR, and mitochondrial imbalance.

A recent study suggested that inhibition of the mitochondrial translation and depletion of respiratory chain subunits would result in the accumulation of large lipid droplets in mouse intestine72. In planarian tissues, phagocytes in the intestine are the major cell type responsible for lipid storage73,74,75. Lipid staining on planarian cryosections also showed a significant increase in lipid droplets in the opa1;egfp RNAi and drp1;egfp RNAi animals compared with opa1;drp1 RNAi and egfp RNAi controls (Fig. 8C, D). To exclude the possibility that larger planarians contain bigger lipid droplets, we compared the worm size and found no significant increase in length after opa1;egfp RNAi and drp1;egfp RNAi compared with opa1;drp1 RNAi and egfp;egfp RNAi controls (Fig. 8E). We then quantified cellular metabolomics data from the four RNAi conditions at 1 dpa (Supplementary Data 9). The results showed that the levels of certain amino acids, including L-citrulline, citrulline, and guanine, decreased in both opa1;egfp RNAi and drp1;egfp RNAi planarians, with a slight recovery in opa1;drp1 RNAi, suggesting a reduction in anabolism when the mitochondrial dynamics was abolished (Fig. 8F). Similarly, metabolites involved in energy metabolism, such as ADP and adenosine 3’,5’-diphosphate disodium salt, increased in both opa1;egfp RNAi and drp1;egfp RNAi planarians and restored in opa1;drp1 RNAi (Fig. 8F), consistent with the reduction of ATP production (Fig. 8A, B). These findings suggest the impaired fatty acid oxidation, in line with the accumulation of lipid droplets in both opa1;egfp RNAi and drp1;egfp RNAi planarians. However, five metabolites, including L-leucine, L-norleucine, glutaryl-carnitine lithium salt, DL-2-aminoadipic acid, and L-2-aminoadipic acid, were significantly increased only in opa1;egfp RNAi and exhibited a recovery in opa1;drp1 RNAi (Fig. 8F). This suggested a unique change in the metabolism after opa1 knockdown. Although a direct target has not been identified thus far, our results suggest the possibility of unraveling the regulation of metabolism by the mitochondrial fusion-fission equilibrium during planarian regeneration.

In summary, we propose a working model that the mitochondrial dynamics can determine the pluripotency of planarian stem cells and is critical for planarian regeneration (Fig. 8G, H). Maintaining the mitochondrial fusion-fission equilibrium is required for the normal transcriptional program in the mitonuclear balance and metabolism (Fig. 8I, J).