Gap Junctional Blockade Stochastically Induces Different Species-Specific Head Anatomies in Genetically Wild-Type Girardia dorotocephala Flatworms

Interestingly, the centroids of wild-type S. mediterranea, D. japonica, and G. dorotocephala morphology groups (the three most closely related species) all vary along a linear axis through the CV space. This reveals the morphological consequences of variation along a canonical variate. Moving from negative to positive values of canonical variate one correlates with the elongation of the auricles, if they exist (with SMs being at the negative end of the spectrum, and GDs and PFs being at the positive end). Moving from negative to positive values of canonical variate two correlates with movement of the auricles, if they exist, from a very anterior position, to a more posterior one (with PFs being at the negative end of the spectrum, and GDs being at the negative end).

Importantly, the patterning changes induced by electric synapse blockade affect the shape of the brain (Figure 4), not only the overall head shape. Worms with pseudo heads that resemble a given species possess brains that resemble that same species, revealing that the functions of this regulatory pathway extends beyond setting the external outline of the head but instead controls the shape of the brain within. These data also suggest a coupling of brain morphology to head shape. Whether brain morphology drives head shape, or whether both are parallel consequences of the same GJC-dependent upstream control mechanism, this finding demonstrates the depth of re-organization of body plan after perturbation of GJC connectivity. Coupling or involvement of the central nervous system in morphogenesis would not be entirely surprising, given past work demonstrating the necessity of long-range neural inputs for correct polarity patterning during regeneration in planaria [20] and vertebrates [109,110].

In addition, perturbation of gap junction communication reorganizes the distribution of the planarian adult stem cells (neoblasts) to spatial patterns appropriate to other species (Figure 5). Together with the change in brain shape, it is clear that the morphological alteration into pseudo morphologies of other species are thorough, involving not only external head shape but also internal organ scaling and positioning of key cell groups. These data reveal the coupling between proliferative stem cells and large-scale morphology, its variability throughout phylogeny, and its control by GJ-dependent signaling. It is possible that control of neoblast distribution is what drives large-scale morphology. Conversely, their spatial distribution may be one of several downstream consequences of the dynamics of a GJC-mediated cellular network.

Our results include the first analysis of the effects of GJ modulation on the long-term bioelectric profile of patterning tissues in vivo. While large-scale bioelectric gradients have been shown to functionally instruct pattern formation in a number of models [12], including determination of morphollaxis/scaling in S. mediterranea [19] and anterior-posterior polarity in D. japonica [19], no studies have yet examined patterns of Vmem distribution in the context of changes in the electrical connectivity of tissues. We found that brief exposure to a GJ blocker, which is known to be washed out of worm tissues within 24 h [20], results in a greater degree of bioelectrical isolation (ability to maintain numerous smaller regions of isopotential Vmem) as late as 10 days after exposure (Figure 6). Although these long-lasting alterations in bioelectrical patterning are not directly caused by octanol remaining within the body of the worm (since 8-OH is known to depart planarian tissues within 24 h [20]), they are initiated by the octanol-mediated disruption of GJC, which alters bioelectric network topology, and results in persistent perturbation of bioelectric state. In this, the effect resembles the known properties of electrical synapses (GJs) in the nervous system, which help implement a kind of Hebbian plasticity that stably alters ionic conductivity among cells after transient electrical modulation [44,111,112,113,114]. Quantitative physiological models of bioelectrical dynamics will need to be developed in order to understand bioelectric state memory [115] and determine why different types of worms have such distinct profiles of regions of isopotential Vmem.

The variability and overall magnitude of the effect imply that the result of 8-OH exposure is only partially effective in shutting down GJC (affecting some percentage of GJs in any animal); this is consistent with the stochasticity of the morphological phenotype (Figure 3). While 8-OH is a commonly-used method for blocking GJC (and in D. japonica, produces the same effect as RNAi targeting innexins genes [20]), we cannot conclusively rule out additional targets besides GJs also being perturbed either directly or secondarily following 8-OH exposure. One difficulty with a molecular (gene-specific) approach to GJ shutdown, is that not only would ~13 innexins need to be cloned out of the GD genome (none have yet been characterized), but their RNAi targeting would have to be tested combinatorially. The space of all possible combinations of different innexin-RNAi’s that would have to be tried simultaneously is enormous, because GJ proteins can often compensate for one another. Nevertheless, our past work tested some such combinations (in D. japonica [20]), and showed that a simultaneous knockdown of three Innexins by RNAi reproduces the 8-OH phenotype, strengthening the conclusion that phenotypes are indeed resulting from 8-OH’s effect on gap junctions.

A recent study [20] showing that transient gap junctional perturbation induces permanent changes in planarian target morphology raised the question of whether this pattern memory was established by alterations of GJC that remain long after the blocking reagent leaves the tissues, or by GJC states that alter biochemical or transcriptional cell states which are stable after GJC returns to normal. While we cannot rule out additional epigenetic changes that might occur, our data show that even 10 days later, the effects of brief 8-OH exposure are preserved as decreased electrical connectivity in the somatic tissues. The ability of GJC networks to stably perpetuate induced changes in coupling are consistent with the known roles of GJs in memory and synaptic plasticity [44,52,114,116,117,118].

Our analysis did not reveal a simple correspondence between the specific distribution of voltage gradients and morphological outcomes of pseudo worm regeneration, which could have several explanations. First, it is possible that not voltage but a different molecular signal is what traverses the GJ network to determine planarian head shape. GJ-permeant molecules with signaling roles include serotonin [119,120] and other neurotransmitters [121], histamine [122], calcium [123], and others, which will be tested in future studies, especially as fluorescent probes for some of these physiological signals come on-line [124,125,126]. Another possibility is that each morphological state is achieved by multiple possible distributions of voltage; in this case, the apparent variability would mask our ability to detect a specific mapping from distribution to head pattern outcome. The most likely explanation however is that shape outcomes may not be encoded via steady (single-image) distributions but in time-dependent activity, as occurs during information encoding in the central nervous system. In this case, movies of voltage change across time in vivo will have to be acquired and analyzed—a significant challenge in the planaria model system. Our analysis of isopotential regions is only a first attempt to mine the information encoded in bioelectric states. Overall, we believe it is very likely that a much more sophisticated analysis method will have to be developed, perhaps akin to what is being used to infer content of neural networks from brain scans [80,127], and applied to a large sample of worm imaging data, to truly crack the bioelectric code in this system.

Functional studies in planaria and Xenopus [19,25,128] have demonstrated the relevance of bioelectrical gradients for anterior patterning in general, but new computational approaches will need to be undertaken to demonstrate what spatio-temporal pattern of bioelectric activity is required for each of the patterning outcomes we observed. Uncovering this mapping will be a focus of our efforts for forthcoming studies.

The pseudo morphologies we observed are not permanent: they eventually remodel, after the first regeneration process is complete, to once again resemble G. dorotocephala. In this, they are distinct from the permanent double-headed phenotype induced in D. japonica, which represents a stable alteration of the worms’ target morphology following physiological perturbation [20,129]. Thus, in G. dorotocephala, repair occurs in two phases: head regeneration following amputation, which completes in the normal timeframe (<10 days) and results in different species’ shapes, and a subsequent much longer remodeling back to the genome-specific G. dorotocephala head shape. This latter process occurs in the absence of injury, within a complete head, over several weeks, illustrating the ability of cellular networks to detect deviations from correct shape, even in the absence of trauma or wounding, and to effect remodeling. Interestingly, the same biphasic process occurs in vertebrates: when a salamander tail blastema is transplanted onto a flank, a tail forms at first, but later remodels into a limb in a much slower process [130,131].

Morphology appears to be a homeostatic parameter, one that is consistently assessed and adjusted in order to maintain some “target” morphology. Morphology must first be decided within days after amputation, as it is possible to observe head shape as early as day six post-amputation. In 8-OH treated worms, this mechanism produces incorrect, pseudo morphologies. However, morphology must then be reassessed again after this first decision. When the current morphology does not agree with the target morphology, remodeling is initiated sometime between day 10 and day 17. The activity of these remodeling processes after regeneration has been completed is fascinating, and implies that morphology is consistently reassessed and edited, even after large-scale re-organization events have ended. The nature of the processes that drive cell behavior towards a specific, stable end-state (when remodeling ceases) is almost entirely obscure. At this time, it is not possible in any model system to derive specific shape information (to know precisely which anatomical pattern will be a sufficient end-goal state to cause remodeling to stop) from genomic data.

The different head morphologies we observed, which are distinct in outward anatomy, distribution of the stem cells, and shape of the brain, must result from differential cell migration and cell-cell communication. To link cell-level behaviors with the discrete large-scale morphologies and explore the consequences of reduced GJ communication (induced by octanol exposure) on emergent cells signaling and morphogenesis, we analyzed an agent-based model. In this model (Figure S2), some cells actively migrate, and alter the shape of the resulting tissues. Their movements are guided by a chemical gradient produced by other static cells in the head, which pattern of activation is established through GJC, as is known to occur in the worm [20,70] and in many other organisms [51,132,133]. Our model shows the morphogenetic consequences of disruption of this communication, which parallels the interplay between signaling of planarian neoblasts (stem cells) and surrounding somatic tissues. In this model, we propose that significant instructive (GJ-dependent) signaling occurs from the somatic tissues to the neoblasts or their progeny, to guide the new tissue generation and shaping during regeneration.

The computational simulations revealed how a small number of cell configurations could result in head shapes quite closely corresponding to the observed morphologies (Figure 9, and Table 2). Interestingly, when specific static cells (black) lose the ability to send their signal due to GJC blockage by octanol, the migrating cells’ movements (red) are altered, giving rise to one of the 4 outcomes. Thus, the model exhibits the two main properties of the experimental dataset: multiple possible outcomes (due to different patterns of GJC-blocked static cells, possibly from physiological noise in pharmacological drug-target interaction), and discrete specific outcomes (resulting from the deterministic signaling of static “black” cells and consequent behavior of migratory “red” cells). It should be noted that our model is also readily applied to an organ that grows from an initial cell position to its final position: if we allow division of cells B, cell migration proceeds in the same way as shown: they will leave behind them deactivated (differentiated) cells B that do not move and do not divide. Future work will test the applicability of this class of models to other regenerating systems.

Our agent-based model includes physical parameters of cells [95,134], as well as traditional chemical gradient diffusion [135,136]. In the model, octanol disruption of cell:cell communication is modeled by the deactivation of signaling from a specific cell type; one candidate for this unique cell population are the rare but widely distributed cells that express the gap junction protein Innexin 5 [71]. Future experiments will test the specific predictions of this model and assign molecular-genetic identities by the several different cell types postulated in the model. An additional important line of investigation is to mechanistically link large-scale dynamics of physiological signaling to the specific behavior of cells modeled in our simulations.

The outcome of regeneration in this model is a set of coherent, discrete head shapes—not a confusing amalgam of different patterns that would be expected if cells made purely local decisions and were merely dis-coordinated by octanol treatment. Currently, the field largely lacks conceptual models for thinking about these kinds of inputs to development and regeneration; it is not clear how best to represent or model discrete, yet stochastically-chosen, emergent large-scale patterns. Yet, this is an increasingly important goal as advances in regenerative medicine and synthetic bioengineering require the ability to predict and control global outcomes by tweaking cell-level parameters. The field of developmental bioelectricity is still searching for an appropriate formalism, but one obvious candidate is that of the neural network paradigm used by computational neuroscientists to understand emergent dynamics of decision-making in brains composed of electrically-coupled cells. We have previously argued that because somatic cells have many of the same components as nerves (ion channels, neurotransmitters, GJs functioning as electrical synapses, gene expression regulated by electrical activity), we may profitably use some of the concepts from neuroscience to understand similar phenomena in slower developmental networks, where output is not muscle-driven behavior but pattern formation [1,13,137].

Scientists in computational neuroscience have modeled the states of electrically-coupled networks in an energy space that reveals the landscape of trajectories that the networks undertake [138,139,140,141,142,143]. Such (neural) networks can stabilize in attractors in this virtual space, and their activity can be parsimoniously explained by networks’ efforts to minimize a global “energy” function [144,145,146,147,148]. Thus, we propose that one way to conceptualize the linkage between bioelectrical connectivity among worm cells (GJC) and discrete global behavioral outcomes achieved by regeneration is to view the worms’ cells as a neural-like network, processing developmental information more slowly but obeying the same statistical dynamics rules as are known to occur in the brain (Figure 9A). It is likely no accident that GJs are well known to be important for information processing in the brain [54,149,150].

We hypothesize that distinct shape outcomes may be attractors in a state space describing distinct modes of the bioelectric network of planarian body cells connected by gap junctions (electrical synapses); in essence, these attractors correspond to regions of a morphospace as classically understood by evolutionary theory [151]. In the experiments, octanol disrupts the normal dynamics of this network, destabilizing and knocking its stochastic dynamics into one of several nearby attractors, normally belonging to different species of flatworm (Figure 9B). Developing a quantitative, computational model of cell networks that have such modes (with electrical states impinging on cell proliferation, apoptosis, and migration behaviors that result in different head shapes) forms the next deep challenge in this area.

Such dynamic networks underlie virtually every complex process known in the biological world, from metabolism to group behavior [152,153,154]. These richly interconnected nonlinear networks, made up of many simple units, have the ability to produce complex pattern in open and thermodynamically unstable systems, without any single governing body [155,156]. Within a given system, these networks converge onto one or more discrete, stable states [157]. Neural networks are just one example of complex dynamic networks; they are implemented via electrical signaling and gap junctions that demonstrate discrete attractors corresponding to global states [141,158], and stable re-occupation of large-scale states (memory) after perturbation [54,147]. Given that neural signaling evolved by specializing much more ancient somatic physiological communication mechanisms, it is reasonable to postulate fundamentally similar network dynamics for physiological networks during pattern formation [137,159]. On this view (Figure 9), the discrete head patterns are outcomes of bioelectric network dynamics with the lowest free-energy; cutting the network raises the free energy by deviating the system from its optimal morphology (damage), and initiates the regenerative process—a journey back down to a low-energy attractor. Perturbing the physiological bioelectric network with 8-OH results in shifting the system dynamics from this attractor to another attractor, and subsequent regeneration of an altered morphology. Apparently, these pseudo shapes are only local minima, since after about 3–4 weeks the pseudo worms return to the global minimum of the standard G. dorotocephala morphology.

In this class of models, the output of the physiological network consists of signals directing cells to proliferate, apoptose, or migrate, so that head morphogenesis is achieved. The details of this process at the cellular level is beginning to be understood, both in the CNS (control of gene expression by brain activity [160]) and developmental bioelectricity [12,23,32]. Variations in membrane potential between cell groups may act as electrophoresing forces, pulling small signaling molecules, small RNAs, or charged transcription factors to specific regions of the planarian body, producing numerous possible long-term effects [47,161,162,163,164]. Additional mechanisms include voltage-sensitive phosphatases [165], and voltage-regulated action of chromatin modifying enzymes [34,43,166]. Clearly, this is just one way to think about these data, and others are possible. The advantage of this paradigm is the rich set of mathematical results available from studies of computational neuroscience that can constrain future model-building and help make testable predictions.

One of the key aspects of these phenotypes is their stochastic nature: the same treatment performed on a cohort of clonal (isogenic) animals results in one of several patterns within any exposed population. Strikingly, the frequencies with which species-specific patterns are observed match the evolutionary distances estimated for these species of planaria (Figure 5). There are significantly more pseudo G. dorotocephala, D. japonica, and S. mediterranea morphologies, which correspond to a smaller evolutionary distance from G. dorotocephala, while P. felina pseudo morphologies remain relatively infrequent. One possibility is that a core network (composed of gene regulatory elements and physiological gap junction communication (GJC)—mediated signals) is at work among many planarians, with species-specific differences tweaking its key dynamics. The evolutionary distance among the planarian types may thus be expected to reflect the ease of shifting this network into its different potential outcome states (as occurs during GJC inhibition). It should be noted that wild-type P. felina have multiple eyes along the doso-ventral body boundary, a trait that our pseudo PFs do not display. Incomplete mimicry of the most distant species, P. felina, is consistent with the association between evolutionary distance and ease with which species-specific morphology can be fully recapitulated during perturbed regeneration. It should be noted that Phagocata gracilis and Phagocata velata planarians also have head morphologies quite similar to P. felina, and therefore also to pseudo PFs. These two species possess only two eyes, and both species are similarly evolutionarily removed from G. dorotocephala. However, no living wild-type worms of these two species could be procured for this work, and therefore our analysis focused on P. felina.

Overall, the ability to generate distinct morphologies from the same genome likely represents an extremely flexible mechanism that could be exploited by evolution via many known environmental and genetically-encoded modulators of gap junctional state.

An obvious absence of “intermediate” morphologies, or head shapes that resemble a blending of two or more species of planarian, could be the result of constraints that limit the morphological space available to the regenerating planarian. Within this framework, the constraints on phenotypic space become driving factors in morphogenesis. These constraints are likely an amalgam of developmental, environmental, and evolutionary factors, which serve to limit the stable modes of the patterning network, and thus the morphogenetic possibilities available to a regenerating flatworm. Efforts currently on-going in our lab are aimed at producing mechanistic bioelectrical models of single cells and cell networks to quantitatively explain the discrete stable modes of such networks, the stochastic nature of each “run” of such a network following amputation, and the molecular links to downstream cell behaviors [94].

Our study faced several limitations. First, we have not yet produced cell-resolution spatio-temporal datasets on voltage and GJ connectivity during remodeling, of the kind used to extract semantics of bioelectrical states in brain activity [80,127]. Second, while electric current is a very likely messenger transferred through the gap junctional paths [19,167], involvement of other small molecules such as neurotransmitters, IP3, and calcium are also possible mediators of the signaling that regulates cell behavior during head morphogenesis [35,119,120,122,168,169]. All of these questions can be investigated once again-of-function technology (misexpression) becomes available in the planarian model system. In addition, future experiments that could reveal differences among pseudo morphologies would target non-coding small RNAs and chromatin modifications.

It should be noted that any model seeking to explain these data would have to include not only a list of molecular components required for the process, but a constructive scheme exhibiting discrete, stochastic, large-scale morphogenetic outcomes. Thus, the biggest area for future development, in addition to molecular details, is computational modeling of a gap junction-based network with sufficient quantitative detail to reproduce the observed stochastic, discrete attractor modes corresponding to specific morphologies. This effort is currently on-going in our lab [137,170,171]. A truly successful model should explain why these particular shapes are produced, and allow the experimenter to control the patterning outcome at will by appropriate manipulation of physiological state and connectivity—a kind of guided self-assembly, as is becoming possible in the frog model system [1,13,17,33].

G. dorotocephala planaria were obtained from Carolina Biological Supply Company. Worms were kept in commercial natural spring water (Poland Springs water, Framingham, PA, USA) at 20 °C, and starved for >7 days before amputation. Pre-tail (PT) fragments were obtained by making transverse cuts posterior to the pharynx and anterior to the tail, using a sharp scalpel. Worms were immobilized by placement on a moistened and cooled Kimwipe. Immediately after amputation, PT fragments were transferred into octanol (8-OH) solution (Sigma Aldrich, RM00050 (St. Louis, MO, USA)), and kept in sterile, deep-dish petri dishes (Fisherbrand, 100 mm × 20 mm). Octanol solution was prepared by slowly pipetting 8 µL 8-OH into 500 mL Poland Springs water, resulting in a final concentration of 123 µM, and allowing the solution to vortex for >30 min. Fragments were then left to regenerate in 8-OH for three days, with daily refreshment of the drug. After three days, fragments were transferred into fresh water until regeneration was complete enough to score head morphology (~10 days).

Worms used for morphometric analysis were imaged using a Nikon SMZ1500 microscope (Melville, NY, USA) with a Retiga 2000R camera (Surrey, BC, Canada) and Q-Capture imaging software (Surrey, BC, Canada). Care was taken to image worms when their heads were most relaxed, in order to capture morphology as accurately as possible (all worms were imaged prior to fixation, which can distort the heads’ external edges). Landmark data were then recorded using ImageJ (Bethesda, MD, USA) [172]. Landmarks for morphometric analysis were chosen based on biological relevance, and reproducibility across worms with varying head morphologies. Landmarks are: (1) the anterior tip of the head; (2) and (3) along the lateral edges at a distance of 1/6th of the length of the worm from the anterior tip of the worm; (4) and (5) the plane of the eyes; (6) and (7) the widest point of the head (often the auricles, if present); (8) and (9) halfway between the eyes and the anterior tip of the head; and (10) and (11) the thinnest point of the head posterior to landmarks 6 and 7. n = 10 WT (wild-type) G. dorotocephala, n = 9 WT D. japonica, n = 6 WT P. felina, n = 8 WT S. mediterranea, n = 7 pseudo G. dorotocephala, n = 13 pseudo D. japonica, n = 5 pseudo P. felina, and n = 6 pseudo S. mediterranea flatworms were measured. For measurements of brain in immunohistochemistry, the posterior edge of the brain was consistently defined as the place where nervous tissue enlarged relative to the ventral nerve cord. MorphoJ software (Manchester, UK) [173] was used for Canonical Variate Analysis, in order to quantify and graphically represent changes in head morphology. MorphoJ software was also used to calculate Procrustes distances, and perform statistical analyses.

DiBAC4(3) (bis-[1,3-dibarbituric acid]-trimethine oxanol) (Invitrogen, Carlsbad, CA, USA) was used as previously [79,83]. A stock solution (1.9 mM) was diluted 1:1000 (0.19 uM) in Poland Springs water, and worms were soaked in the DiBAC solution for >30 min before imaging. Worms were then immobilized in 2% low-meting point agarose, using custom-fabricated Planarian Immobilization Chips as in [174]. Images of the ventral side of immobilized planaria were captured with the Nikon AZ100 Stereomicroscope, Melville, NY, USA, using epifluorescence optics, and NIS-Elements imaging software (Melville, NY, USA). Data were neither altered nor removed.

Whole worm immunohistochemistry was done as in [175]. Primary antibodies were 3C11 anti-synapsin raised against mouse (used at 1:50 dilution), obtained from the Developmental Studies Hybridoma Bank, created by the NICHD of the NIH and maintained at the University of Iowa, Department of Biology, and α-phosphorylated histone H3 (H3P) 1:250 (Upstate). Secondary antibody was goat anti-mouse Alexa488 (Sigma, St. Louis, MO, USA), used at 1:400 dilution (for anti-synapsin), and an HRP-conjugated anti-rabbit antibody with TSA-Alexa568 and anti-HRP (Molecular Probes, Eugene, OR, USA) (for anti-H3P).

In order to estimate the degree of gap junctional connectivity, we processed voltage dye images to determine the number of isopotential cell fields in MATLAB (MathWorks Inc., Natick, MA, USA). Images were flat-field corrected in Image J (Bethesda, MD, USA). For each image, the area of the flatworm was found by taking the average intensity of all of the pixels in the entire image. The flatworm pixel area was used to calculate five intensity clusters of isopotential for the worm body. Each pixel of the flatworm was categorized into the closest of the five clusters by Euclidean distance. For each image, all of the pixels categorized into one intensity cluster were isolated by thresholding, regions of size less than 15 pixels were removed, and any image holes were filled. The boundaries of the region were then calculated using the Moore-Neighborhood tracing algorithm with a modified Jacob’s stopping criteria (MATLAB function bwboundaries) [176]. The regions from each intensity cluster for each flatworm were summed to compute the total number of regions in an image (Figure 6I).

Numerical simulations were carried out with an agent-based model where cells are considered as mass points interacting with each other. Some cells have fixed positions and some cells can move. The forces acting between cells depend on the distance between them. Cell motion is calculated as the motion of their centers according to Newton’s second law. The forces acting on cells include friction/adhesion and pressure/chemotaxis from other cells. Cells can be added, moved, and activated/deactivated by the user to create the initial configuration and during a simulation, so any initial cell configuration can be established. The membrane configuration is fixed during a simulation. Figure S2 contains additional details about the modeling. Numerical modeling was carried out with original C++ code developed for this work. The new software we developed includes a user interface with parameter windows and a computational domain (Figure S3).

Morphological remodeling poses immense computational challenges for the multicellular organism. Cellular state must be integrated and communicated across time and anatomical distance; parameters of regenerative or developmental programs must be set, and these programs must be implemented to narrow tolerances. While the genome is the source of all structural material in the organism, gap junctions provide unique conduits for using physiological and bioelectrical signals to link cells into information-processing networks, and are therefore ideal complements to genetically-encoded programs. We show that a relatively simple physiological perturbation can derive distinct species-specific morphologies from the same genome, highlighting an important layer of control between the genetics and the anatomy of an organism. Further work is needed to better characterize both the dynamics of these networks and the molecular mechanisms that they inform, as well as develop tools with which to control morphology in experimental contexts and biomedical applications.