Genome-wide analysis reveals conserved transcriptional responses downstream of resting potential change in Xenopus embryos, axolotl regeneration, and human mesenchymal cell differentiation

Endogenous bioelectric signaling via changes in cellular resting potential (V mem) is a key regulator of patterning during regeneration and embryogenesis in numerous model systems. Depolarization of V mem has been functionally implicated in dedifferentiation, tumorigenesis, anatomical re‐specification, and appendage regeneration. However, no unbiased analyses have been performed to understand genome‐wide transcriptional responses to V mem change in vivo. Moreover, it is unknown which genes or gene networks represent conserved targets of bioelectrical signaling across different patterning contexts and species. Here, we use microarray analysis to comparatively analyze transcriptional responses to V mem depolarization. We compare the response of the transcriptome during embryogenesis (Xenopus development), regeneration (axolotl regeneration), and stem cell differentiation (human mesenchymal stem cells in culture) to identify common networks across model species that are associated with depolarization. Both subnetwork enrichment and PANTHER analyses identified a number of key genetic modules as targets of V mem change, and also revealed important (well‐conserved) commonalities in bioelectric signal transduction, despite highly diverse experimental contexts and species. Depolarization regulates specific transcriptional networks across all three germ layers (ectoderm, mesoderm, and endoderm) such as cell differentiation and apoptosis, and this information will be used for developing mechanistic models of bioelectric regulation of patterning. Moreover, our analysis reveals that V mem change regulates transcripts related to important disease pathways such as cancer and neurodegeneration, which may represent novel targets for emerging electroceutical therapies.

Keywords: Axolotl, depolarization, differentiation, embryogenesis, ion channel, mesenchymal stem cells, microarray, transcriptome, Vmem, Xenopus

Along with biochemical signals, cell behaviors such as migration, morphological change, proliferation, and differentiation are regulated by physical properties. In particular, bioelectric signaling among non‐neural cells has recently been shown to be an instructive component of pattern formation during regeneration, development, and cancer (Borgens 1988; McCaig et al. 2005; Zhao et al. 2006; Forrester et al. 2007; Stewart et al. 2007; Zuberi et al. 2008; Chernet & Levin 2013a; Yang & Brackenbury 2013). One important aspect of endogenous bioelectricity is resting potential or V mem (Levin 2014b). V mem controls proliferation, differentiation, apoptosis, and migration of a wide range of cell types (Gilbert & Knox 1997; Wang et al. 1999; Weihua et al. 2005; Blackiston et al. 2009; Sundelacruz et al. 2009), including mammalian stem cells (Sundelacruz et al. 2008, 2013a; Pillozzi & Becchetti 2012; Swayne & Wicki‐Stordeur 2012; Shen et al. 2013; Sundelacruz et al. 2013b; Wang et al. 2014). Moreover, spatio‐temporal patterns of V mem have been shown to specifically regulate growth and form (Levin 2012, 2013, 2014a; Adams & Levin 2013; Tseng & Levin 2013a). Ion channel activity, and the resulting voltage gradients, regulate shape, size, and positional information of organs in Drosophila, planaria, fish, frog, salamander, and mouse (Adams et al. 2007; Ozkucur et al. 2010; Lange et al. 2011; Dahal et al. 2012; Beane et al. 2013; Perathoner et al. 2014), as well as being responsible for several classes of developmental malformation in humans (Galanopoulou 2010; Tristani‐Firouzi & Etheridge 2010; Masotti et al. 2015). Bioelectric signals can trigger the formation of whole ectopic organs, such as in the case of eye development (Pai et al. 2012a), and can stimulate the repair of complex tissues such as spinal cord (Borgens et al. 1990, 1999; Shapiro et al. 2005; Tseng et al.

2010). Endogenous bioelectric gradients can also control the morphogenesis of limbs, faces, and whole body axes (Tristani‐Firouzi & Etheridge 2010; Marrus et al. 2011; Vandenberg et al. 2011), making these gradients an important modality for targeted intervention in regeneration and bioengineering applications (Levin & Stevenson 2012).

Despite the extensive functional data summarized above, there remains limited insight into specific pathways and mechanisms that directly link V mem to such phenotypic changes. Recent work has identified several transduction mechanisms by which voltage potential change is converted into chromatin modification and changes of expression of a handful of specific target genes (Levin 2012; Tseng & Levin 2012) in several contexts. However, no systematic analysis has been performed to identify and integrate genome‐wide transcriptional changes following steady‐state depolarization in vivo, or to suggest novel biomedical endpoints for V mem modulation. Moreover, while bioelectricity has been investigated in many phyla, the evolutionary conservation of downstream signaling pathways has not been fully described. For example, does the bioelectric code (the mapping of voltage states to specific organ outcomes) function via the same transcriptional programs in diverse organism morphologies, or is there a significant divergence in the downstream responses to bioelectrical signaling? Answering such questions is crucial, not only to fully understand the developmental role of physical forces, but also to help establish roadmaps for human regenerative medicine from data on bioelectric control of patterning in model systems of development and repair. Such roadmaps would also provide a path forward for the systematic probing of the pathways and the control achievable via V mem, analogous to what has been pursued for biochemical signaling mechanisms over the past few decades.

Here, we use genome‐wide analysis of transcriptomes as a tool for gaining insight into transcriptional cascades, and the degree to which these are conserved, downstream of V mem change. Such strategies have been applied in several model organisms (Altmann et al. 2001; Baldessari et al. 2005; Chalmers et al. 2005; Tomancak et al. 2007; Yanai et al. 2011) and assist in acquiring a holistic view with systematic characterization of how gene expression is regulated towards functional paths during processes such as development or reproduction (Martyniuk & Denslow 2012; Langlois & Martyniuk 2013). In this study, we analyzed the effects of specific depolarization events on transcriptional profiles in Xenopus laevis embryos during development. We then compared these data with those from the regeneration of spinal cord in axolotl, Ambystoma mexicanum, and human mesenchymal stem cells (hMSCs). The strategy was to identify common features of signaling downstream of V mem change along two orthogonal dimensions: species (frog, salamander, and human), and cellular context (embryogenesis, regeneration, and stem cell differentiation in vitro).

Regulation of resting potential (V mem) of a wide variety of cell types plays an important role in embryogenesis, regenerative response, and cancer (Borgens et al. 1977b; McCaig et al. 2005; Blackiston et al. 2009; Pullar 2011; Lobikin et al. 2012; Adams & Levin 2013; Chernet & Levin 2013b, 2014; Tseng & Levin 2013b; Levin 2014a; Pai & Levin 2014). In order to understand the role of bioelectricity in pattern formation, and to harness this signaling modality for biomedicine, it is important to understand the transcriptional networks downstream of specific V mem change. The microarray experiments reported here (from the three different species Xenopus laevis, axolotl and human) represent single time point experiments, performed to enable direct comparison of multiple transcriptomics datasets. In the following discussion we focus on only some of the biological themes affected by depolarization; our decision was guided by our research interests and paper‐length limitations. However, all the data for processes affected by depolarization are provided in the Appendices, to enable others to pursue any other interesting leads.

A recent study of transcriptional changes during Xenopus embryonic development identified dynamic changes in membrane hyperpolarization gene networks which were decreased in early development but increased later in development at stage 34 (Langlois & Martyniuk 2013). Here, we used microarray analysis to identify transcripts that are regulated specifically by depolarization (induced by the activity of each of two very different depolarizing channels, as is done during bioelectric induction of patterning changes in vivo). We focused on a single time point, just after mid‐gastrula transition, when new transcription begins (Woodland & Gurdon 1968; Forbes et al. 1983; Cascio & Gurdon 1987). Further detailed studies with different time points during development will be explored in future work, to understand long‐term, temporal aspects of V mem regulation of gene regulatory networks during embryonic development.

Changes in endogenous V mem patterns encode important signals for individual cells and also for large‐scale patterning programs. To begin to analyze the global transcriptional targets of such change, we induced specific depolarization in frog embryo cells by misexpression of a depolarizing ion channel. We previously showed that this technique results in a specific, coherent anatomical change (induction of well‐formed ectopic eyes throughout the animal [Pai et al. 2012a]), and it is a strategy routinely used to investigate bioelectric signaling (Adams & Levin 2013; Adams et al. 2013). Xenopus laevis embryos were injected with either 666 (DN‐KATP) (Hough et al. 2000), glycine‐gated chloride channel (GlyR) (Davies et al. 2003), or water (controls). The mRNA extracted from each of these treatments (n = 50 each) was used for microarray analysis using an Affymetrix X. laevis Genome Genechip 2.0 Array (Fig. 1). DN‐KATP has been previously shown to cause depolarization of the injected cells in Xenopus embryos by inhibiting KATP channels (Hough et al. 2000; Pai et al. 2012a). Similarly, expression of the GlyR channel in the presence of the channel opener drug ivermectin (IVM) also depolarizes the injected cells in Xenopus embryos (Davies et al. 2003; Blackiston et al. 2011; Pai et al. 2012a). We used two different (K+ and Cl− ion flux) channels that both depolarize embryonic cells, in order to focus on genes whose transcription is specifically responsive to depolarization, not sodium or potassium signaling per se (nor on any possible ion‐independent functions of one channel protein).

In the Xenopus dataset, 380 genes were significantly upregulated in both experimental groups (DN‐KATP and GlyR+IVM), in comparison to controls (water‐injected), with at least a twofold increase in gene expression (Appendix S1). There were 140 genes that were commonly downregulated in both experimental groups in comparison to controls, with at least a twofold decrease in gene expression (Appendix S1). Since the list of downregulated genes was small, we focused our analysis on the significantly upregulated gene list from the Xenopus dataset. To organize our data into functional categories, we used subnetwork enrichment analysis (SNEA). These networks were then verified with the results of a functional classification analysis using the PANTHER database.

To begin to understand the degree of conservation of these transcriptional responses, we also performed a microarray analysis of axolotl spinal cord regeneration. Control (vehicle water) treated animals were compared with animals treated with the depolarizer IVM (opens endogenous GlyR channels leading to Cl− efflux). In these experiments, the IVM or control vehicle (water) was directly injected into the central canal of the axolotl spinal cord, and the spinal cord was injured by removing a portion of the spinal cord (Diaz Quiroz & Echeverri 2012; Sabin et al. in review) (Appendix S1). Tissue samples were harvested 1 day post‐injury, and tissue from 10 axolotls was pooled and RNA extracted for microarrays. Arrays were carried out using a commercially available custom Affymetrix axolotl array. In the axolotl dataset, 756 genes were significantly upregulated in the experimental groups (IVM treatment 1 day post‐injury), in comparison to controls (water‐injected), with at least a twofold increase in gene expression (Appendix S1). There were 753 genes that were commonly downregulated in the experimental group in comparison to controls, with at least a twofold decrease in gene expression (Appendix S1). We did not restrict the bioinformatics analysis by fold change, and have included all genes that showed an uncorrected P < 0.05 (Appendix S1).

To extend the analysis to human stem cells, as well as to compare the above in vivo assays with results obtained in vitro, we compared the above microarray data to microarray analysis of osteoblasts derived from undifferentiated hMSCs. Normal osteoblast differentiation medium has low extracellular K+. The osteoblasts were depolarized by (1) elevating extracellular K+ via addition of potassium gluconate (40 mmol/L) into the medium (elevating the extracellular K+ reverses the electrochemical gradient for K+ leading to depolarization of cells) or (2) treatment with the Na+/K+ ATPase inhibitor ouabain (OB) (10 nmol/L) in the medium. Depolarization induced by these concentrations of OB and K+ has been confirmed using voltage‐sensitive dyes and/or sharp intracellular recordings (Sundelacruz et al. 2008; and Kaplan DL et al., unpublished data). Osteoblasts in the normal differentiating medium were used as controls (Sundelacruz et al. 2013a). The microarray was performed using Illumina Human WG6 v3 Expression BeadChip arrays, which have 48,804 probe sets, of which over 27,000 represent coding transcripts with well‐established annotation. We used two different treatment conditions that both depolarize cells, in order to focus on genes whose transcription is specifically responsive to depolarized groups, n = 3, and a non‐depolarized osteogenic control group, n = 3, were used. In the human dataset, 2777 genes were significantly upregulated in both the experimental groups (high K+ in medium and OB treatment), in comparison to controls (normal differentiation medium) (Appendix S1). There were 2706 genes that were downregulated in comparison to controls (Appendix S1). We however did not restrict the bioinformatics analysis by fold change, and have included all genes that showed a P < 0.01 (Appendix S1).

In summary, we reasoned that common pathways identified in the comparison (across species and across process type) represent those cell signaling pathways and processes most probably affected following membrane depolarization. However, additional efforts must continue to verify the role of changing V mem on these pathways.

SNEA was performed for each of the three datasets, an analysis that determines if there are significantly enriched processes in a list (based on gene function) compared to a background list (annotated subnetworks in the database based on mammalian literature) (Figs. 2−6, S1−S3, Tables 1, 2, S1 and S2). Common subnetworks were identified, and this was followed by manual grouping into major biological themes.

Out of a total 3441 cell processes in the database (Fig. 2A, List 4), there were 100, 284, and 372 cell processes affected in Xenopus, axolotl, and human cells respectively following depolarization (Fig. 2A). There were 2981 cell processes unaffected by any treatments across all three experiments (Fig. 2A, red). Similarly, out of a total 4545 disease networks in the database (Fig. 2B, List 4), there were 83, 469, and 113 diseases affected in Xenopus, axolotl, and human cells respectively following depolarization (Fig. 2B). There were 4065 disease networks unaffected by any treatments across all three experiments (Fig. 2B, red). It is interesting to note that a total of 69 cell processes and 51 disease networks were enriched in common between the frog, axolotl and human datasets in response to depolarization (Fig. 2A, B and Tables 1 and 3). Each of these networks contained a minimum of 10 entities/genes. Comparing these to the total list of possible cell processes and disease networks (Fig. 2A, B) revealed that only a subset of cell processes and disease networks were affected by depolarization at the level of the transcript.

Of the common processes, there was approximately 25% that were related to organ development (17/69) while 23% were related to cell cycle and cell function (16/69). All subnetworks are presented in Appendix S2. Some examples of interesting depolarization‐affected gene networks involved in organ development are depicted in Figure 3. Noteworthy is that the gene networks regulated by depolarization include organ/organ systems belonging to all three embryonic germ layers (ectodermal, mesodermal, and endodermal) (Figs. 3, S1 and Table S1). This was found to be the case in the other two (axolotl and human) datasets as well (Table 1). Transcripts involved in brain/neural development (Fig. 3B), skeletal/bone development (Fig. 3C), immune system development (Fig. S1), and adipose development (Fig. S1) were found to be significantly affected by depolarization.

V mem is known to regulate cell behaviors, although the transcriptional mediators of voltage‐regulated shape change are poorly understood (Ghiani et al. 1999; Bauer & Schwarz 2001; Chittajallu et al. 2002; Sundelacruz et al. 2008; Blackiston et al. 2009; Becchetti & Arcangeli 2010; Schwab et al. 2012; Stock et al. 2013). To address this, we searched for “large‐scale functions” (regulated processes that are found across all tissues and organ structure and are necessary for proper developmental morphogenesis and anatomical sculpting) which may be regulated by depolarization. We identified several large‐scale function gene clusters within the frog dataset (Fig. 4B, Table 1). Depolarization‐mediated regulation of such large‐scale functions was also found to be a common theme across all three (frog, axolotl, and human) datasets (Table 1). Our analysis confirmed V mem regulation of key cell cycle events in vivo during development (G2/M transition, mitotic entry, S‐phase and DNA replication) (Table 1 and Fig. 4). Our analysis identified putative transcriptional targets by which V mem regulates apoptosis (Table 1 and Fig. 4), and also identified other cell death mechanisms like anoikis that are regulated by V mem (Table 1). Our analysis also identified, for the first time, fate specification genes regulated by V mem for tissues from all three germ layers (Table 1) (e.g., nkx2.5, gata4, nkx3.1, hnf1b) (Tronche & Yaniv 1992; Shiojima et al. 1995; Bhatia‐Gaur et al. 1999; Zhu et al. 2000; Kohler et al. 2008; Li et al. 2012; Zhou et al. 2012).

We were interested in identifying cellular‐level signaling pathways that are responsive to depolarization, particularly those that are conserved across species. SNEA was performed for “expression targets” in Pathway Studio. Indeed, a common subset was found to be differentially affected in all three datasets (Table S2). During embryonic development a collection of conserved juxtacrine and paracrine signals (consisting of growth factors, morphogens, hormones, and cytokines) are known to integrate together in various spatial and temporal permutations−combinations to drive specific developmental patterns collectively defined here for the first time as a “developmental simaton”. The majority of the common subset of pathways identified in our experiments belonged to the developmental simaton (Table S2). These include IGF, FGF, BMP/TGFβ, HGF, EGF, PDGF, and gonadotropins among others. SNEA for “expression targets” revealed that BMP2 may be an important developmental simaton pathway regulated by depolarization (Figs. 5A, S2A, C and Appendix S3) due to a significant number of transcripts with altered expression in the datasets that were downstream of this signaling molecule in all three species. Depolarization also regulates cell cycle (oncogenic/tumor suppressor) signals across species (Table S2). Another set of depolarization‐regulated signals conserved across species (frog, axolotl, and human datasets) are ion translocators, particularly those mediating calcium and chloride transport and signaling (Figs. 5B, S2B, C). Some examples of depolarization‐regulated ion channel (particularly calcium and chloride) genes in the Xenopus dataset are voltage‐dependent calcium channel p/q type alpha 1A (cacna1a), voltage‐dependent calcium channel beta unit 4 (cacnb4) and pendrin (slc26a4).

To determine which diseases might be associated with depolarization, we used SNEA to query disease networks associated with differential gene expression in all three (frog, axolotl, and human) datasets. The largest group of diseases were neoplasms associated with various tissues (Table 2). Surprisingly, a significant number of metabolic disease networks like diabetes, insulin resistance, and glucose intolerance were also regulated by depolarization (Fig. 6A and Table 2) as were neural disorders (Figs. 6B and S3).

In order to confirm the observations seen by SNEA, we performed an independent analysis of the frog dataset using the PANTHER (protein annotation through evolutionary relationships) functional classification system (Mi & Thomas 2009; Mi et al. 2013a) (Fig. 4A and Table 3). This algorithm is designed to classify genes and their proteins to specific functions. We used the “biological process” function, which clusters the genes and their functions in the context of larger networks. Since a gene can be classified according to more than one term, the pie chart is calculated according to the number of “hits” to the terms divided by the total number of “class hits.” A class hit indicates an independent (not parent or child to each other) ontology term (Mi & Thomas 2009; Mi et al. 2013a). Developmental process was one of the 13 biological processes enriched from the entire gene list of upregulated genes. PANTHER analysis also suggested that genes upregulated by depolarization in the frog dataset were involved in regulating embryonic developmental processes that spanned across all three germ layers (Fig. 4A), similar to the SNEA observations (Tables 1, S1 and Fig. 4). These processes included anatomical structure morphogenesis, pattern specification process, ectoderm development, endoderm development, mesoderm development, and systems development (Fig. 4A). A deeper analysis of genes clustered within the “systems development” category showed that the regulated genes were involved in the development of organ/organ systems, for example muscle development, nervous system development, skeletal system development, and heart development (Fig. 4A). These categories are in line with those obtained using SNEA (Tables 1, S1 and Fig. 3).

The PANTHER analysis supports the observation that depolarization‐mediated changes in mRNA abundance affect organogenesis of tissues and organs across all three germ layers during embryonic development.

We independently analyzed the pathways in the frog dataset using the PANTHER database classification system (Mi & Thomas 2009; Mi et al. 2013a) (Table 3). Developmental simaton pathways and cell cycle regulatory (oncogenic/tumor suppressor) pathways were a major subset of the pathways regulated by depolarization similar to SNEA (Tables 3 and S2). Neural, neurological disorders and immune pathways were some other pathways that overlapped with the SNEA results (Tables 2, 3, Figs. 6, S1 and S3).

We investigated the evolutionary conservation of transcriptional responses to depolarization across two dimensions: species (comparing between frog, axolotl, and human), and process type (development, spinal cord regeneration, and stem cell differentiation). Our analysis included single time points, not a temporal time course. Comparative analysis for the effect of V mem signals was done under similar V mem conditions (depolarization). Thus, we reasoned that common pathways identified in the comparison represent those cell signaling pathways and processes most likely affected following membrane depolarization. However, additional efforts must continue to verify the role of changing V mem on these pathways. Below we discuss a subset of genes and processes regulated by V mem (selected by us based on the likely impact in several key fields of biology); a complete list is provided in Appendices S1−S4. We focused on Xenopus because it is in this system that the largest number of published mechanistic data are available, demonstrating that endogenous V mem gradients and their modulation can be used to achieve predictable and coherent changes in body pattern formation. Hypotheses derived from our analysis will be most readily first tested in frog, but then moved into human tissues and other model systems. As such, the discussion remains focused on this species with reference to the other two models as a comparison.

To identify cell processes and diseases related to genes affected by depolarization events, differentially expressed genes from all three datasets (Xenopus, axolotl, and human) were queried in Pathway Studios using SNEA. Comparing these to the total list of possible cell processes and disease networks (Fig. 2A, B) revealed that only a subset of cell processes (460 out of 3441, ∼13%) and disease networks (480 out of 4545, ∼10%) were affected by depolarization at the level of the transcript. The above data on large‐scale global changes in gene expression pattern reveal that the effect of membrane potential depolarization is not indiscriminate or random. While additional functional experiments are required, transcriptomics suggest that depolarization may affect specific cell processes and diseases (either directly or indirectly, via secondary genetic or biochemical signals), commonly across the three species. A total of 69 cell processes and 51 disease networks were enriched in common between the frog, axolotl, and human datasets in response to depolarization (Fig. 2A, B and Tables 1, 2), suggesting that these are the downstream events most probably affected by depolarization. The major biological theme underlying depolarization‐mediated changes in cell processes appears to be the regulation of tissue and organ development during embryogenesis (17/69, 25%, related to organ development, and 16/69, 23%, related to cell function and cell cycle) (Table 1 and Appendix S2). This includes depolarization‐mediated regulation of cellular level processes (cell fate, cell death, cell interactions), tissue level processes (tissue size, function, and interaction), organ level processes (development of muscle, neural, skeletal, and heart/vascular tissues), and organ‐system level processes (embryonic development processes, organ architectures, regeneration) (Table 1).

It is important to note that the number of cell processes affected is higher in the human dataset but the number of disease processes affected is higher in the axolotl dataset. This is probably due to the differences in the transcriptomes of specific cell types in culture (mesenchymal) versus an entire spinal cord region including the surrounding cells and tissues involved in regeneration (with contributions from all three germ layers). The disease pathways for the human dataset derived from a single germ layer and homogenized cell culture are limited compared to the axolotl dataset derived from all germ layers and numerous kinds of cells in that tissue. It is a conservative approach to analysis, as a much larger list of commonalities would be expected if a comparison were to be made with a human tissue sample similar to the axolotl sample with contributions from all three germ layers.

Organogenesis consists of complex series of spatial and temporal patterning events. Since our data analyzed transcriptional state at only one point, it is revealing only a snapshot of the responses to depolarization. However, previous studies analyzing Xenopus developmental processes at multiple time points show V mem as a critical regulator of organogenesis at several stages (Langlois & Martyniuk 2013). To further understand what aspects of organogenesis are affected by depolarization, we analyzed the frog dataset for regulation of organ development processes using SNEA. Noteworthy is that the gene networks regulated by depolarization include organ/organ systems belonging to all three embryonic germ layers (Fig. 3 and Table S1). This was found to be the case in the other two (axolotl and human) datasets as well (Table 1).

V mem coordinates differentiation and morphogenesis of ectodermal craniofacial (Vandenberg et al. 2011) and eye and brain structures in X. laevis embryogenesis (Pai et al. 2012a, 2012b, 2015). SNEA identified brain/neural tissue patterning as significantly regulated processes by V mem, supporting the previous observations. However, this analysis suggests that specific processes within brain/neural tissue patterning are regulated by V mem signals (synaptic transmission, neuroprotection, neurite outgrowth), and identifies the genes regulated by V mem in each of these processes (Fig. 3B and Table 1). Similar V mem‐mediated regulation of target genes involved in neural/brain tissue development and disease was also observed in axolotl spinal injury and human cells (Figs. 3B, 6B, S3 and Table 1). This analysis generates testable hypotheses for discovering the precise role and mechanism of V mem in regulation of neural tissue patterning, especially with respect to bioelectrical control of the NCAM, BMP, FGF, serotonergic, adrenergic, and dopaminergic pathways.

In addition to ectodermal tissues (brain, eye, and craniofacial structures), this analysis for the first time suggests that V mem also regulates mesodermal (heart, muscle, bone, adipose, immune system) and endodermal (kidney, pancreas, lung) organogenesis (Fig. 3 and Tables 1 and S1). Gene networks related to skeletal/bone development were identified (Fig. 3C) based on the differentially expressed transcripts. These included osteoclast differentiation, osteoblast differentiation, ossification, biomineralization, calcium mobilization, and calcium‐mediated signaling. Analogous to neural/brain development, the gene targets significantly enriched in skeletal/bone development subprocesses and the adipocyte differentiation process (Fig. S1) are found to overlap with similar gene‐target‐enriched subprocesses regulated by depolarization in axolotl and human datasets (Table 1). These observations are further supported by recent studies showing V mem‐mediated regulation of adipogenic and osteogenic differentiation of hMSCs (Sundelacruz et al. 2008, 2013a), and the induction of regeneration of bony appendages induced by bioelectric signaling (Borgens et al. 1977a, 1977b, 1979, 1983, 1984; Sisken et al. 1984; Tseng & Levin 2013a).

Endoderm mainly forms the tubular lining of internal organs associated with digestive tract (gut, liver, pancreas, and kidneys) and lungs (Fukamachi & Takayama 1980; Kedinger et al. 1990; Gilbert 2010). We found evidence of V mem regulation of endodermal genes and processes: in Xenopus, the hnf1b gene is involved in kidney and pancreas development and atmin is involved in kidney and lung tubulogenesis (Table S1). Also V mem regulation of endodermal processes such as gastrulation, tubulogenesis, lumen formation, branching morphogenesis, and hepatic regeneration (Gilbert 2010) was found to be a common significantly enriched gene target theme in all three organisms (Table 1). The limited observed effect of V mem on transcripts involved in endoderm organogenesis could be a consequence of the temporal limitations of this study, since endoderm development programs were not highly active at the time point of sample collection.

PANTHER analysis of the frog dataset also showed that “developmental processes” was one of the 13 enriched biological processes and those subsets of genes were involved in regulating embryonic developmental processes that spanned across all three germ layers (Fig. 4A), similar to the SNEA observations (Tables 1, S1 and Fig. 4). Moreover, the genes were involved in muscle, skeletal, nervous, and cardiac/heart system development (Fig. 4A), categories in line with those obtained using SNEA (Tables 1, S1 and Fig. 3). The PANTHER and SNEA together indicate that depolarization‐mediated transcriptomic changes affect organogenesis of tissues and organs across all three germ layers during embryonic development.