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

Tissues and organ systems develop through regulation of basic large‐scale functions such as proliferation (size control), cell differentiation (fate determination), dedifferentiation (cancer and regeneration), and cell death (morphogenetic sculpting) (Abud 2004; Stanger 2008a, 2008b). These functions collectively direct changes in both growth and form to generate the appropriate anatomical pattern. Functional data have shown that bioelectric state change could reverse or duplicate whole body axes, including the head−tail axis in planaria (Oviedo et al. 2010; Beane et al. 2011) and the left−right axis in frog, chick, and zebrafish embryogenesis (Levin et al. 2002; Adams et al. 2006). How is cell membrane depolarization able to alter large‐scale anatomy? We identified several large‐scale‐function gene clusters within the frog dataset (Fig. 4B, Table 1). These included regulation of various stages of the cell cycle, cell survival, cell proliferation, cell differentiation, dedifferentiation, and cell death. 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). These data suggest that depolarization regulates development of tissues and organ/organ systems across all germ layers, perhaps by controlling such basic large‐scale functions found across all cells and tissues.

The rhythmic oscillation of V mem throughout the cell cycle (Bregestovski et al. 1992) occurs in concert with the regulation of cell cycle and proliferation by stable resting potential in a range of embryonic, somatic, and neoplastic cells (Blackiston et al. 2009; McCaig et al. 2009; Sundelacruz et al. 2009; Lobikin et al. 2012; Adams & Levin 2013; Chernet & Levin 2013b, 2014). In this paper we identified cell cycle genes like fos, p27, and mapk, which were previously shown to be regulated by V mem (Kong et al. 1991; Wang et al. 2003), and we also identified novel links between V mem and some classic cell cycle regulatory genes such as tp53 and Ras. These last two are particularly interesting because they suggest the presence of feedback loops in the control of oncogene‐driven tumorigenesis by depolarization (Chernet & Levin 2013a, 2013b, 2014).

In addition to being potent drivers of cell differentiation along various lineages (forming several different kinds of cells) (Ghiani et al. 1999; Bauer & Schwarz 2001; Chittajallu et al. 2002; Sundelacruz et al. 2008), V mem signals also possess the capacity to induce dedifferentiation and transdifferentiation (Cone & Tongier 1971; Harrington & Becker 1973; Stillwell et al. 1973; Cone & Cone 1976; Sundelacruz et al. 2008, 2013a). In vivo experiments have shown the remarkable power of V mem signals to re‐specify organ identity by inducing cell differentiation across germ layers (Beane et al. 2011, 2013; Pai et al. 2012a, 2015). Here we identified for the first time fate‐specification genes regulated by V mem for tissues from all three germ layers (Table 1 and Appendix S2) (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). These data suggest that other organs such as heart, liver, kidney, pancreas, bone fat, muscles, and immune system may also be inducible by appropriate V mem modulation, in addition to the ectopic brains, eyes, and limbs that have been produced so far.

Programmed cell death or apoptosis mediates tissue morphogenesis during development and is in fact required for regeneration (Schwartz 1991; Abud 2004; Tseng et al. 2007; Li et al. 2010; Ryoo & Bergmann 2012). Recent studies have shown V mem regulation of apoptosis as a way of regulating morphogenesis and tissue remodeling (Lang et al. 2005; Beane et al. 2013; Englund et al. 2014). But how V mem regulates apoptosis is not well understood. Our analysis not only identified putative transcriptional targets by which V mem regulates apoptosis (Table 1 and Fig. 4), but also identified other cell death mechanisms such as anoikis that are regulated by V mem (Table 1). It offers valuable insight into mechanisms of cell death regulation by indicating mitochondrial membrane potential and genes such as bax as targets of V mem signals (Table 1 and Appendix S2). The harnessing of these endpoints by targeted bioelectric modulation, using currently available technologies such as pharmacological cocktails (Tseng et al. 2010; Famm et al. 2013; Sinha 2013), optogenetics (Bernstein et al. 2012; Adams et al. 2013; Spencer Adams et al. 2014), and gene therapy (Chernet & Levin 2013b, 2014; Pai et al. 2015), is an exciting area for the extension of synthetic bioengineering techniques.

A common subset of cellular‐level signaling pathways responsive to depolarization was found to be differentially affected in all three datasets (Table S2). The majority of these pathways belonged to the “developmental simaton” (defined here for the first time as a conserved collection of juxtacrine and paracrine signals for growth factors, morphogens, hormones, and cytokines known to integrate together in various spatial and temporal permutations−combinations to drive specific developmental patterns) (Table S2). SNEA for “expression targets” in Pathway Studio showed that BMP2 was regulated in the frog, axolotl, and human cell experiments (pre‐differentiation) (Figs. 5A, S2A and Appendix S3). BMP2 was also identified as being enriched as a regulator of downstream expression targets that was also differentially regulated in the dataset (i.e., gene hub). Our findings support recent work showing BMP expression to be downstream of ion channel function (Dahal et al. 2012; Swapna & Borodinsky 2012). Strikingly, a deeper analysis of the Xenopus array revealed that the developmental simaton (conserved across species) pathways such as insulin‐like growth factor receptor 1 (igfr1), fibroblast growth factor receptor 3 (fgfr3), and others are involved in depolarization‐regulated tissue/organ development processes across all three germ layers (brain/neural, skeletal/bone, adipose, and immune system), supporting their role as a master‐regulator for specific organs/organ systems (Tseng & Levin 2013a; Levin 2014a). Here we have analyzed only a single snapshot of the interaction of V mem with the developmental simaton. A temporal analysis is required for revealing further dynamic interactions between V mem and components of the developmental simaton.

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). Previous studies have shown that calcium signaling is an important transduction mechanism for V mem‐mediated regulation of brain and eye patterning (Beane et al. 2011; Pai et al. 2012a). Interestingly, there appears to be crosstalk between the developmental simaton signals like BMP2 and BDNF and the ion‐translocator‐mediated signals (Figs. 5B, S2B, C). The regulation of calcium and chloride transport and signaling, and their crosstalk with the developmental simaton signals, suggests the presence of feedback loops, where the bioelectric state change regulates (directly or indirectly via its influence on the developmental simaton) expression of new ion channels which in turn further alter the bioelectric landscape an ongoing cycle of feedback between physiological and transcriptional dynamics.

Another mechanism downstream of voltage change is bioelectric regulation of serotonergic signaling, which has been characterized in the role of V mem in left−right patterning, neural pathfinding, and melanoma‐like transformation (Levin et al. 2006; Blackiston et al. 2011, 2015; Lobikin et al. 2012). Our analysis indicated additional neurotransmitters which have not yet been investigated in this context (Table 3 and Fig. S3). One example is slc1a3 the sodium‐dependent glutamate/aspartate transporter which suggests that these neurotransmitters should be tested in developmental bioelectricity assays. Similar observations have also been made with respect to the dopamine system in Xenopus development (Langlois & Martyniuk 2013).

Overall, the data indicate that there are certain signaling pathways induced by depolarization that seem to regulate processes at all levels of organization cellular, tissue, organ, and organ system levels. Specifically, depolarization regulates cell cycle pathways (oncogenic/tumor suppressor) which control large‐scale (basic broad spectrum) functions like cell fate, dedifferentiation, proliferation, and death. It can be hypothesized that V mem control of such basic large‐scale function, in turn, may allow it to regulate development of multiple tissues and organs across all germ layers during embryonic development. This is supported by the documented role of bioelectricity as a master‐regulator for specific organs/organ systems (Tseng & Levin 2013a; Levin 2014a). Our data suggest that bioelectric regulation of developmental simatons may be one of the mechanisms by which ionic signaling exerts long‐range, non‐cell‐autonomous effects not only during developmental patterning (Vandenberg et al. 2011; Pai et al. 2012a, 2012b, 2015) but also in tumor suppression (Chernet & Levin 2014), metastatic induction (Blackiston et al. 2011; Lobikin et al. 2012), and regenerative remodeling (Oviedo et al. 2010; Beane et al. 2011, 2012, 2013; Lobo et al. 2012). Most importantly, the various targets of depolarization are found to be largely conserved themes among three very diverse model systems (frog, axolotl, and human) and distinct in vivo and in vitro contexts (development, regeneration, and stem cell biology) analyzed. This suggests a conserved set of responses that will facilitate the medical exploitation of bioelectric control methods.

Based on gene expression profiles, our data revealed that oncogene/tumor suppressor genes are the second largest group of pathways regulated by depolarization (Tables 3 and S2), providing a mechanistic link between the known involvement of ion channels in cancer (Kunzelmann 2005; Pardo et al. 2005; Felipe et al. 2006; Stuhmer et al. 2006; Prevarskaya et al. 2010; Lobikin et al. 2012; Yang & Brackenbury 2013; Lang & Stournaras 2014). To determine which other diseases might be associated with depolarization, we identified the disease networks associated with differential gene expression in all three (frog, axolotl, and human) datasets. As expected, the largest group of diseases were neoplasms associated with various tissues (Table 2). Numerous studies have now shown an intimate relationship between ion channels and human cancers, and targeting ion channels is now being considered as a novel therapy for personalized non‐genetic cancer treatment (Schonherr 2005; Fiske et al. 2006; Arcangeli et al. 2009; House et al. 2010; Lobikin et al. 2012).

Surprisingly, a significant number of metabolic disease networks like diabetes, insulin resistance, and glucose intolerance were also related to genes regulated by depolarization (Fig. 6A and Table 2). Other disease networks include those related to neural diseases, immune disorders, cardiac disorders, pulmonary disorders, developmental disorders/birth defects, and wound healing and regenerative defects (Table 2). We observed an underlying theme of hypertrophic diseases across tissues in relation to depolarization, for example cardiac hypertrophy, cardiomegaly, pulmonary fibrosis, pulmonary hypertension, and even inflammation and arthritis (Table 3). This is probably due to the regulation of basic cell processes like cell proliferation, differentiation, and death by depolarization. While a range of channelopathies reveal ion channels at the root of birth defects (see Table 1 in Levin [2013]), most of these have not been specifically linked to depolarization. Thus, our data suggest specific clinical endpoints that should also be investigated with respect to resting potential change. Our long‐term goal is to apply this knowledge to develop bioelectric strategies to address disease states for regenerative medicine.

The neural disorders, especially neurotoxicity and nerve degeneration (Figs. 6B and S3), are particularly interesting as we found neurotransmitter pathways (Table 3) and neurodegenerative pathways (e.g., Parkinson's and Huntington's) (Table 3 and Fig. S3) to be induced by depolarization, suggesting that resting potential changes may represent a functional therapeutic target in these cases. Experiments in Xenopus have shown that enforcing correct nervous system V mem signals can rectify developmental brain defects (Pai et al. 2012a, 2015). Such information is valuable because a large panel of ion channel drugs exists, many already approved for human use, and can be tapped as off‐label uses for electroceuticals in the nervous system and non‐neural tissues. Another disease group of note is the wound healing and regeneration defects, providing gene targets as hypotheses to be tested during V mem’s known effects on wound healing and regenerative processes (Chifflet et al. 2005; Cao et al. 2011a, 2011b; Vieira et al. 2011; Luxardi et al. 2014).

A large number of cell processes that are regulated by V mem change were found to be in common among amphibian and human systems, and across contexts of embryogenesis, spinal cord regeneration, and adult stem cells (Fig. 7). Our data are consistent with a highly conserved role of bioelectric state as a major regulator of cell activity and pattern regulation at all levels of organization cellular, tissue, organ, and whole body axes, both in vitro and in vivo. The consequences of depolarization were observed in transcripts characteristic of all germ layers (ectoderm, mesoderm, endoderm, and neural crest), highlighting resting potential as a powerful control point for biomedical interventions. The identification of neurotransmitter pathways linked to voltage change in somatic cells suggests not only new models of signal transduction but also the possible relevance of models of synaptic information processing that may be relevant to cellular decision‐making during pattern regulation. The identification of disease pathways suggests that bioelectric signaling and intervention using electroceutical strategies should be examined for diseases such as cancers of several different tissues, metabolic disorders like diabetes, congenital malformations (e.g., neural tube defects), and neurodegenerative diseases. Taken together, our data provide a framework for testing mechanistic hypotheses for bioelectric pathways, as novel components for synthetic bioengineering circuits using ionic signaling, and as targets for biomedical intervention in a range of disease states.

Xenopus laevis embryos were fertilized in vitro according to standard protocols (Sive et al. 2000) in 0.1× Marc's modified Ringer's solution (10 mmol/L Na+, 0.2 mmol/L K+, 10.5 mmol/L Cl−, 0.2 mmol/L Ca2+, pH 7.8). Intracellular ion concentrations in Xenopus embryos are 21 mmol/L Na+, 90 mmol/L K+, 60 mmol/L Cl−, 0.5 mmol/L Ca2+ (Gillespie 1983). Xenopus embryos were housed at 14–18°C and staged according to Nieuwkoop and Faber (1967). All experiments were approved by the Tufts University Animal Research Committee (M2014‐79) in accordance with the guide for care and use of laboratory animals.

Capped synthetic mRNAs generated using mMessage mMachine kit (Ambion, ThermoFisher Scientific, Grand Island, NY, USA) were dissolved in nuclease‐free water and injected into the embryos in 3% Ficoll using standard methods (Sive et al. 2000). Each injection delivered between 1 and 2 nL or 1 and 2 ng of mRNA (per blastomere) into the embryos, usually at the one‐cell stage into the middle of the cell in the animal pole. Constructs used were GlyR (Davies et al. 2003) and 666 chimera (Hough et al. 2000).

The procedure was performed by the Beth Israel Deaconess Medical Center (BIDMC) Genomics Core (Boston, MA) at Harvard University. A miRNeasy Mini kit (Qiagen, Hilden, Germany) was used to isolate total RNA from the Xenopus embryos. There were 50 embryos pooled for each experimental sample. Microarray hybridization target preparation was performed using the Affy 3′IVT Express Kit (Affymetrix, Santa Clara, CA) as per the manufacturer's protocol. Fragmented and biotin labeled/amplified RNA was hybridized to the GeneChip® Xenopus laevis Genome 2.0 array (Affymetrix, Santa Clara, CA) as per the protocol provided by the manufacturer. The Affymetrix GeneChip® X. laevis Genome 2.0 Array has 32,400 probe sets representing more than 29,900 X. laevis transcripts. The quality of hybridized arrays was assessed using Affymetrix guidelines on the basis of scaling factor, background value, mean intensity of chip and 3′ to 5′ ratios for spike‐in control transcripts. The outlier analysis was performed using unsupervised clustering and principal components analysis. All high quality arrays were normalized using the MAS5 algorithm developed by Affymetrix. The absent/present calls for the transcripts were calculated using the MAS5 algorithm. The differentially expressed transcripts were identified on the basis of fold change and Affymetrix transcript calls. Overexpressed transcripts were at least twofold in the experimental group compared to the control group with present call in the experimental group. Under‐expressed transcripts were at least twofold changed in the experimental group compared to the control group with present call in the control group. Differentially expressed genes were annotated using the Affymetrix database or by performing BLAST analysis. All microarray data were analyzed using Bioconductor packages in R. The NCBI GEO accession number for frog data is GSE72099.

Each (RNA) sample, including the control, was generated not from a single individual embryo but by pooling 50 embryos (including control). This gives a robust average expression compared to using single embryos, and smooths out any individual variability. In addition, we pursued a more nuanced and strict replication strategy than simple replicates; this was necessary here because of the specific nature of depolarization as a developmental signal. Traditional replicates would filter out technical noise but would still include genes specifically responsive to a given experimental technique (e.g., IVM exposure), as opposed to what we wanted to study: the effects of depolarization per se. We used two different depolarization methods/samples (each one with 50 embryos pooled together), and kept only those transcripts that were shown to be similarly regulated in both treatments. This new strategy not only filters out technical noise (as would traditional replicates), but also filters out transcripts that are not consistent between the two different ways of depolarization a much stronger criterion that removes a lot more noise, both technical variability and responses that are specific to a single type of perturbation. This provides a stringently curated set of targets related to what we wanted to study depolarization and comes from two methods/samples (each sample representing 50 embryos).

All axolotls used in these experiments were bred in the axolotl facility at the University of Minnesota under the IACUC protocol #1201A08381. Axolotls of 2–3 cm were used for all in vivo experiments, and animals were kept in separate containers and fed daily with artemia; water was changed daily. Animals were anesthetized in 0.01% p‐amino benzocaine (Sigma‐Aldrich, St. Louis, MO, USA) before microinjection was performed.

IVM or vehicle only (water) was pressure injected into the central canal of the spinal cord, and this was visualized by the addition of Fast Green into the solution. Directly after injection, a portion of the spinal cord was surgically removed and the animals were placed back into water in individual containers. One day post‐injury animals were anesthetized again and the area of the injury was removed. Tissue from 10 animals was pooled for each microarray replicate.

Total RNA extraction for microarray analysis was done using TRIzol® Reagent. RNA was resuspended in 20 μL of RNAase‐free water. RNA concentration was measured using a Nano‐Drop 2000 spectrophotometer (ThermoFisher Scientific, Grand Island, NY, USA). RNA integrity was evaluated using an Agilent 2100 Bioanalyser. Microarrays were carried out in triplicate using a custom‐made axolotl Affymetrix Chip, and arrays were processed at the microarray facility at the MPI‐CBG, Dresden, Germany. Briefly, for each sample, 150 ng total RNA was used to reverse transcribe double‐stranded cDNA and subsequently in vitro transcribe biotin‐labeled target cRNA as per the GeneChip 3′IVT Express Kit. The target cRNAs were hybridized to Amby002: a custom Affymetrix GeneChip©. The chip contains 20,000 axolotl probes. Hybridizations were also performed according to the GeneChip 3′IVT Express Kit as per the user manual. Approximately 12.5 μg fragmented and labeled aRNA was hybridized for 17 h at 45°C and 60 rpm. Arrays were washed and stained using the GeneChip Fluidic station FS450 with the wash and stain protocol FS450_0001. The arrays were scanned using the Affymetrix GCS3000 System (Scanner) using default parameters to obtain background‐corrected signal intensities. Cel files containing normalized intensity data (RMA normalization) were generated using the Gene Expression Console (Affymetrix). FDR was not used for the axolotl data. The NCBI GEO accession number for axolotl data is GSE72099.

Whole bone marrow aspirate from a 25‐year‐old healthy man was purchased from Lonza through their Research Bone Marrow Donor Program, following approved guidelines of informed consent as previously documented (Sundelacruz et al. 2008, 2013a). Aspirate was plated at a density of 10 mL of aspirate per square centimeter in control medium (Dulbecco's modified Eagle's medium with 10% fetal bovine serum, penicillin [100 U/mL], streptomycin [100 mg/mL], and 0.1 mmol/L nonessential amino acids) supplemented with basic fibroblast growth factor (1 ng/mL) (Invitrogen, ThermoFisher Scientific, Grand Island, NY, USA). Cells were maintained in a humidified incubator at 37°C with 5% CO2. The hMSCs were isolated on the basis of their adherence to tissue culture plastic and were used for experiments between passages two and four.

Undifferentiated hMSCs were cultured in control medium. Osteogenic differentiation medium consisted of α‐modified minimum essential medium supplemented with 10% fetal bovine serum, penicillin (100 U/mL), streptomycin (100 mg/mL), 10 mmol/L β‐glycerophosphate, 0.05 mmol/L l‐ascorbic acid‐2‐phosphate, and 100 nmol/L dexamethasone (Sigma‐Aldrich).

V mem was depolarized by (1) addition of OB (10 nmol/L; Sigma‐Aldrich) to differentiation medium or (2) elevation of extracellular K+ by adding potassium gluconate (40 mmol/L; Sigma‐Aldrich) to differentiation 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; Kaplan DL et al. unpublished data). Osteogenic cells were pre‐differentiated for 3 weeks before the addition of OB or potassium gluconate for an additional 3 weeks of culture. The depolarizers were added to the differentiation medium and replenished in subsequent media changes.

RNA was isolated from samples collected in TRIzol® reagent according to the single‐step guanidinium acid‐phenol method. RNA was further purified using the RNeasy Mini kit according to the manufacturer's instructions, which included an on‐column DNase treatment. RNA integrity was evaluated with the Agilent Bioanalyzer. The cDNA reverse transcription, cDNA purification, in vitro transcription of cRNA, and microarray hybridization, staining, and scanning were performed by the Yale Center for Genomic Analysis. We used 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. Sample sizes were as follows: depolarized groups, n = 3; non‐depolarized osteogenic control group, n = 3; undifferentiated control group, n = 2.

Data analysis was performed using Bioconductor software packages. The lumi package (Du et al. 2008; Lin et al. 2008) was used to perform variance stabilizing transformation and quantile normalization, and the limma package (Smyth 2004) was used for linear modeling and differential expression analysis. Log2‐fold gene expression changes and empirical Bayes moderated t statistics were computed, with P values adjusted according to the Benjamini−Hochberg method and considered significant for P < 0.01. We performed hypergeometric testing on Gene Ontology terms using the GOstats package (Gentleman 2004) and adjusted P values for multiple testing, considering results significant for P < 0.01. The NCBI GEO accession number for hMSCs data is GSE72099.

SNEA was performed in Pathway Studio 9.0 (Elsevier Life Science Solutions) and ResNet 9.0 to construct gene interaction networks for transcripts showing differential expression in frog, axolotl, and human cells. For each of the three datasets, the list of differentially expressed genes was mapped into the program using official gene symbols (Name + Alias). SNEA was performed and significantly enriched processes were determined to be those with P < 0.05 that also contained more than 10 members in the network. These were the networks that are most represented by the entities in each gene list. Both cell processes and disease subnetworks were queried and there were 500 permutations of the data to generate the distributions. Briefly, SNEA uses known relationships among genes (e.g., relationships based on co‐expression patterns, binding, or involvement in common pathways) to build networks focused around gene hubs. These interaction maps are generated using information from the ResNet 9 database. The database contains over 20 million PubMed abstracts and ∼2.4 million full‐text articles (22 September 2014). Thus, these are pre‐defined molecular networks based on the literature (i.e., it is the background or reference group). A gene list is then imported into the program and a statistical comparison between the experimental subnetworks mapped to known networks and the entire background of known networks (reference group) is conducted using a Mann−Whitney U test; a P value is generated that indicates the statistical significance of difference between two distributions (additional details on the method can be found in the technical bulletin, page 717, from Pathway Studios 7.0). Venn diagrams (Oliveros 2007) were generated using both official gene symbols and pathways to identify common genes and pathways affected in all three experiments.

Gene networks based upon expression, binding, and regulatory interactions of entities were constructed using direct connections with one neighbor. In order to simplify the pathways, only those interactions with the highest scores (most well supported by the literature) are shown (>3000 connectivity; local connectivity >5). However, it is important to note that the pathways are built with additional information that is not shown in the figures (i.e., those connections showing <3000 connectivity). Subnetworks that included “cell process,” “expression targets,” and “disease” were queried for all three species. We reasoned that those processes in common among the datasets were those most responsive to depolarization. These data were then subjected to manual categorization into “major biological themes.”

The gene lists from the microarray experiments were uploaded to the PANTHER database and functional classification was viewed as pie charts as previously documented (Thomas et al. 2006; Mi et al. 2013a, 2013b).

We thank Asipu Sivaprasadarao for the DN‐KATP plasmid, Daryl Davies and Miriam Fine for the GlyR plasmid, Manoj Bhasin and Marie Bruno‐Joseph for assistance with the microarrays, and Joan Lemire for comments on the manuscript. We gratefully acknowledge the support of NSF (EBICS sub‐award CBET‐0939511), the G. Harold and Leila Y. Mathers Charitable Foundation and the AHA (14IRG18570000), the NIH (RO1 AR005593; R01 AR061988), and the W. M. Keck Foundation. The authors declare no conflict of interest.