Positional plasticity in regenerating Amybstoma mexicanum limbs is associated with cell proliferation and pathways of cellular differentiation.

BACKGROUND: The endogenous ability to dedifferentiate, re-pattern, and re-differentiate adult cells to repair or replace damaged or missing structures is exclusive to only a few tetrapod species. The Mexican axolotl is one example of these species, having the capacity to regenerate multiple adult structures including their limbs by generating a group of progenitor cells, known as the blastema, which acquire pattern and differentiate into the missing tissues. The formation of a limb regenerate is dependent on cells in the connective tissues that retain memory of their original position in the limb, and use this information to generate the pattern of the missing structure. Observations from recent and historic studies suggest that blastema cells vary in their potential to pattern distal structures during the regeneration process; some cells are plastic and can be reprogrammed to obtain new positional information while others are stable. Our previous studies showed that positional information has temporal and spatial components of variation; early bud (EB) and apical late bud (LB) blastema cells are plastic while basal-LB cells are stable. To identify the potential cellular and molecular basis of this variation, we compared these three cell populations using histological and transcriptional approaches.
RESULTS: Histologically, the basal-LB sample showed greater tissue organization than the EB and apical-LB samples. We also observed that cell proliferation was more abundant in EB and apical-LB tissue when compared to basal-LB and mature stump tissue. Lastly, we found that genes associated with cellular differentiation were expressed more highly in the basal-LB samples.
CONCLUSIONS: Our results characterize histological and transcriptional differences between EB and apical-LB tissue compared to basal-LB tissue. Combined with our results from a previous study, we hypothesize that the stability of positional information is associated with tissue organization, cell proliferation, and pathways of cellular differentiation.

Background

Urodele amphibians such as salamanders and newts are exceptional model organisms to study processes of endogenous reprogramming and regeneration because they are capable of regenerating complicated biological structures from mature adult tissues. Regeneration in these organisms occurs through the modification of mature adult cells into regeneration-competent cells, known as blastema cells, which re-pattern and re-differentiate into the missing structures. One fascinating aspect of the regenerative process is that the regenerating blastema tissues only replace the structures that are missing, and thus a unique developmental pattern is established in the blastema depending on the location of the injured tissue. There are two main hypotheses to explain when and how this unique pattern is established in regenerating limbs, both of which have compelling supporting evidence. The first hypothesis is known as the “pre-specification model”, and is based on the idea that the pattern of the entire regenerate is established in the blastema as soon as it forms [1]. The second hypothesis, which has multiple sub-hypotheses, is broadly based on the idea that the pattern of the regenerate is generated gradually as the blastema develops [2, 3].

Recently published observations from our lab and others, suggest that patterning information is initially plastic and is gradually stabilized in specific regions of the blastema as it develops. One piece of evidence that strongly supports this hypothesis is that positional information in blastema cells can be reprogrammed when exposed to Retinoic Acid (RA) [4, 5]. For example, if an early bud (EB) stage blastema is exposed to RA, its positional information is reprogrammed such that it will form a complete limb regardless of the location of the blastema along the proximal/distal axis of the limb [4, 6]. A different reprogramming phenotype occurs if a late bud (LB) stage blastema is exposed to exogenous RA. In this situation, positional information of cells at the apical tip of the blastema is reprogrammed, while cells in the basal region closest to the stump are unaltered [4]. We attribute this difference in reprogramming capacity to reflect the stability of positional information in these tissues. We suspect that the state of tissue differentiation may be related to the stability of positional information in these cells because basal blastema cells are actively differentiating, while apical cells remain undifferentiated. One observation that supports this hypothesis is that positional information in uninjured limb tissues is unaffected by RA exposure [4, 5], despite the expression and activation of RA receptors in these cells [7, 8]. These observations reveal that an EB blastema and the apical tip of a LB blastema are capable of positional reprograming, and that differentiated cells are resistant to being reprogrammed.

While the above-described observations on RA-treated blastemas are consistent with the idea that EB and apical-LB cell populations are positionally plastic, we have further tested this hypothesis by determining whether positional reprogramming could occur through endogenous interactions between these and mature tissues at different positions along the limb [9]. Our experiments rest on the following observation – when cells with stabilized positional information are confronted with cells with differing positional information, the confronted cells induce the formation of new pattern (resulting in the formation of new structure) [10, 11]. Thus, if the confrontation of cells from different locations on the limb does not induce the formation of new structure, then one or both of the populations do not have stabilized positional information. Because mature limb tissues have stabile positional information [9, 12, 13], then the ability to induce the formation of new structure must reflect the stability of positional information in the “other” population of cells. Within this logical framework, we grafted EB and apical-LB tissues next to mature tissues at different locations along the limb. These confrontations did not induce the formation of new structures [9]. Since EB and apical-LB blastema tissues express region-specific Hox genes [1, 2], it is likely that they have positional information. Thus, the inability to induce new structure when grafted to a different position on the limb is probably due to plasticity (or instability) of positional information, rather than lack of this information. Also consistent with this idea is the observation that the expression of region-specific T-Box (Tbx) genes is altered in EB grafts such that they match the expression profile of the new host location [9].

As we stated above, uninjured limb tissue has stabile positional information, and thus can induce the formation of new structures when confronted with cells with differing positional information. We have found that the basal region of the LB blastema also has this inductive capacity, revealing that positional information is stabile in the basal region of a LB blastema. Thus, a LB blastema is composed of both positionally plastic (apical) and positionally stable (basal) populations of cells [9].

Understanding the nature of positional plasticity and the mechanisms that control this property is not only important for understanding pattern formation during limb regeneration, but also could help us improve the efficacy of regenerative therapies that attempt to engraft cells with positional information [14]. With these objectives in mind, we used histological and microarray analysis to identify cellular processes, genes, and signaling pathways that differentiate EB and apical-LB cell populations, from the basal-LB population. Relative to the EB and apical-LB populations, the basal-LB population showed greater actin cytoskeleton and ECM structural organization, and yielded higher expression estimates for genes that promote cellular differentiation. Combined with our previous results [9], we hypothesize that that positional plasticity may be associated with factors that regulate ECM organization and cellular differentiation.

Results

Characterization of tissue organization in EB and apical-LB compared to basal-LB populations

The recent identification of blastema populations with varying degrees of positional stability led us to better characterize these cells and identify differences in cell behavior. We looked at the organization of the actin-cytoskeleton within the cells of the blastema mesenchyme, and the extracellular matrix (ECM) surrounding these cells (Fig. 1). The actin cytoskeleton was analyzed on sagittally oriented tissue sections from EB and LB stage blastemas that had been stained with phalloidin-rhodamine, and the degree of order (alignment) or disorder of actin filaments was quantified using automated image processing that computed the discrete entropy of the actin fibers (Fig. 1b, c). We found that the actin fibers in the populations of cells that are positionally plastic (i.e., EB blastema and apical-LB blastema) were short in length and disorganized such that they did not have any apparent polarity (Fig. 1b, c). The extracellular matrix surrounding these tissues, visualized by immunoflourescence staining of the ECM molecule tenascin, also appeared disorganized (Fig. 1d). In contrast, the average entropy of actin fibers in the basal-LB was lower (Fig. 1b, c), and the extracellular matrix appeared to be more organized (Fig. 1d). We note that similar differences in the density and apparent organization of the reticular lamina, which is composed of multiple ECM molecules including collagens, was observed in the apical and basal regions of limb blastema by Neufeld and Aulthouse [15]. Thus, our observations on tenascin organization likely represent the blastema ECM as a whole. In summary, the amount of tissue organization was relatively low in EB and apical-LB compared to basal-LB populations.

Fig. 1

Characterization of the actin cytoskeleton and extracellular matrix (ECM) during regeneration. a A DIC image of a late bud blastema transverse section indicating where the sagittal sections in (b, d) are located along the proximal/distal axis. b Phalloidin-rhodamine staining of sagittal sections through an EB, or from the Apical to basal regions of the LB. c Quantification of the organization of the actin cytoskeleton in images in (b). Error bars are the SEM, and students t-test was used to determine statistically relevant changes in organization (p < 0.002). d Immunodetection of the ECM molecule tenascin (Green) in EB, or apical to basal regions of the LB

Characterization of cell proliferation among blastema cell populations

According to the Polar Coordinate model of regeneration, when cells with positional information from different locations in the limb are juxtaposed in a regeneration-competent environment, they induce a growth response that generates new cells with the missing positional information, a process known as intercalation [10, 11]. Thus, according to this model, cell proliferation is an integral part of the intercalary response during regeneration. To gain insight on when and where an intercalary response could be occurring in the developing blastema, and how this corresponds to the populations of cells with plastic or stabile positional information, we analyzed cell proliferation in EB blastemas and apical-LB blastemas. The proliferating cells in EB blastemas and LB blastemas were analyzed by calculating the percentage of cells in different regions of the blastema mesenchyme that incorporated EdU into newly synthesized DNA (Fig. 2). We found that the labeling index throughout the mesenchyme of the EB blastema was approximately 20 % (Fig. 2b). At the LB stage, the overall amount of cell proliferation was higher, and the labeling index varied greatly depending on whether the cells were located apically or basally. At the apical-most tip of the LB blastema, approximately 50 % of the cells were labeled. The labeling index increased to above 80 % in the cell population just proximal to these apical-most cells. The labeling index then gradually decreased in the basal region of the blastema (Fig. 2b). Similar results have been reported in regenerating amphibian limb blastemas and developing mouse limb buds [16, 17]. The labeling index in the stump tissue closest to both EB and LB blastemas was around 12 % (Fig. 2b); this is lower than the overall labeling index of the blastema mesenchyme, but higher than the labeling index that has been reported in uninjured limb tissues [18, 19].

Fig. 2

Characterization of cell proliferation in early and late blastemas. a Fluorescent images obtained on 10uM transverse sections of 7 day (early) and 15 day (late) blastemas that had been stained with DAPI (blue) for nuclei and Edu (green) for dividing cells. Animals were injected with Edu 3 h before tissue was harvested. White lines indicate the boundaries of each region that was quantified (100 uM segments for 7D samples, and 200 uM segments for 15D samples). b Histogram of the labeling index of digital sections from the apical (left) to basal (right) region of each type of blastema. The labeling index was quantified in 4 complete blastema sections for each sample. Error bars are the SEM, students t-test was used to determine statistically significant differences (p < 0.05) in the labeling index of the most basal blastema data-points and the stump data-points (red-filled symbols)

Together, these results show that EB and apical-LB blastema populations have increased amounts of cell proliferation relative to the neighboring stump and basal-LB populations, respectively. It is possible that the increased amount of proliferation in EB and apical-LB tissue is associated with the process of intercalation. This idea is explored further in the discussion.

Analysis of EB, apical-LB, and basal-LB blastema transcriptomes

To better understand mechanisms underlying the property of positional information in the regenerate, we tested for transcriptional differences among EB, apical-LB, and basal-LB blastema tissues that developed from proximal (mid-humerus) and distal (carpals) sites of limb amputation. We included proximal and distal blastema tissues in the experimental design because results from our previous study [9] predict that positionally plastic and positional stable cell populations should exhibit similar transcriptional patterns for genes associated with positional information, regardless of the site of amputation. Using a custom axolotl Affymetrix array (20,080 probesets) [20], 365 probesets were identified as significantly differently expressed using a two-way ANOVA model that treated blastema sample type (EB, apical-LB, basal-LB) and amputation site (proximal, distal) as fixed effects (Additional file 1: Table S1). To visualize these results in the form of a heatmap (Fig. 3), we performed two clustering analyses: 1) An objective methodology (see Methods and Additional file 2: Figure S1) was used to group the differently expressed probesets into 15 clusters; 2) Then, differently expressed probesets were used to cluster the blastema samples according to gene expression similarity. The interaction term (amputation site x blastema sample) was significant for 163 of the 365 probesets, indicating that gene expression varied among blastema populations as a function of amputation site. For example, 42 of the significant probesets in Cluster 1 were significant for the interaction term. Typically, Cluster 1 genes were expressed more highly in the proximal basal-LB sample than in all other samples, while expression estimates for proximal apical-LB and EB samples were similar to all three distal samples. As another example of gene x amputation site interaction, genes in Cluster 9 were similarly expressed in proximal basal-LB and apical-LB, but estimates tended to be higher and similar for proximal EB and all three distal samples. These results show that transcript abundances for many genes varied as both a function of apical-basal location within the blastema and proximal-distal location of the amputation site.

Fig. 3

Significantly different gene expression in EB, apical-LB, and basal-LB mesenchymal tissues. The heatmap displays gene clusters (y-axis) and sample clusters (x-axis). Genes and samples were both clustered by average linkage clustering using 1-Pearson’s correlation as a distance measure. Only genes that tested statistically significant for the overall model were used for cluster analysis. Fifteen gene clusters were chosen and colored using R’s rainbow function [72]. Each colored box represents gene expression in a single sample for a single probe set with yellow indicating high expression and red indicating low expression. Samples grouped according to cell population, with the first major division occurring between proximal basal samples (samples B2, B3, B5, and B6) and all other samples. The next major division separates proximal apical samples from the remaining samples, and the following division separates proximal samples from the basal samples. Within distal samples, EB and apical-LB samples clustered together

Only 24 probesets were significant for the amputation site main effect; the majority of these (N = 16) were assigned to Cluster 2. In comparison, 178 probesets were significant for the blastema sample main effect and these were distributed among all clusters except Cluster 12. As was observed for the list of significant interaction probesets, many of these genes were assigned to Clusters 1 and 9. The probesets from Cluster 1 (N = 49) registered transcript abundance estimates (within amputation sites) that were higher for basal-LB samples and lower for EB and apical samples. In contrast, probesets from Cluster 9 registered relatively higher expression estimates for EB and apical-LB samples than the basal-LB sample.

Identifying a prioritized list of significant genes

Based on our previous study [9] and assuming that transcript abundances correlate with properties of positional information, we identified a high priority list of candidate genes (Additional file 3: Figure S3). We used t-tests (p < 0.05) to identify the set of probesets within the overall list of 365 that did not exhibit significantly different transcript abundances between proximal or distal EB and apical-LB tissues (N = 127). We then filtered this list to identify probesets that yielded statistically significant transcript abundances between EB and apical-LB compared to basal-LB tissues. To accomplish this filtering step, we used t-tests (p < 0.05) to identify significant probesets for the following four contrasts: Significantly different expression between EB and basal-LB blastemas arising from (1) distal and/or (2) proximal amputations; significantly different expression between apical-LB and basal-LB blastemas arising from (3) distal and/or (4) proximal amputations. A total of 67 significant probesets were identified for two or more of the contrasts (Additional file 4: Table S2). These results show that genes associated with developmental patterning (bicd2, osr1), chromatin structure (jarid2, nfia), cell proliferation, growth, and differentiation (hbp1, cdkn2c, cndp2, gas2), and major signaling pathways (FGF- fgfrl1, fgfr3, fgfbp3; RA- crabp1; BMP/TGFb- tgfbr2, crim1; WNT- tcf7l2) are expressed differently when comparing EB and apical-LB to basal-LB blastema.

Annotation of significant genes

Of the 365 significant probesets, 325 had assigned gene names. We used these gene names to test for significantly enriched biological process ontology terms using Panther gene analysis tools [21] and a conservative familywise error rate (p < 0.0001). We ran separate analyses for each of the clusters in Fig. 3, however only Clusters 1 and 9, the two largest gene clusters, yielded significantly enriched terms. Eight Cluster 1 genes (acan, sox6, fgfr3, osr1, gdf5, tgfbr2, lect1, myf5), which were more highly expressed in basal-LB samples, significantly enriched the cartilage development ontology term (Bonferoni corrected prob = 0.008). These and other Cluster 1 genes are predicted to pattern tissues during early development, including structures associated with skeletal (e.g., fgfrl1, foxc1, col9a2) and nervous systems (rgmb, rpe65, acan, fat3, sox6, fgfr3, gdf5, col9a2, tgfbr2, crim1, tcf7l2, pcdh15, gfra2, dact1, cdkn2c, gpr37L1, nrcam, ptprd, fabp7, lrig1, rgma, epha5, foxc1) (Additional file 1: Table S1). A different set of related ontology terms were identified for Cluster 9 genes, which were more highly expressed in apical-LB and EB samples. Several ontogeny terms associated with extracellular matrix organization were enriched by an overlapping set of genes, including serpinh1, mmp1, mmp2, mmp13, sparc, ctsk, fn1, itgad, nfkb2, ero1l (Additional file 1: Table S1) (Bonferoni corrected prob = 0.004). These genes encode proteins associated with matrix structure, disassembly, and collagen catabolism.

To further explore the significant gene list, we searched the literature using gene names as queries. We focused on genes involved in cell signaling and chromatin modification because both would seemingly be required to induce and maintain a plastic state. In Table 1, we highlight genes that fall within four general categories: cell signaling, chromatin modification, cell metabolism, and neural function/development. These genes encode proteins that function in FGF, ESRRG (estrogen-related receptor gamma), and mechanotransduction signaling pathways, as well as genes that modify histones via methylation, acetylation, and ubiquination.

Table 1

Genes with higher expression in basal-LB tissues

Potential role in blastema	Abbreviation	Gene name	Reference
Regulation of tissue differentiation	CAPSL	Calcyphosine-like	[77]
	CRABP1	Cellular retinoic acid binding protein 1	[78]
Skeletal differentation and development	ACAN	Aggrecan	[79]
	COL9A2	Collagen, type IX, alpha 2	[80]
	DACT1	Dishevelled-Binding Antagonist Of Beta-Catenin 1	[42]
	EMILIN3	Elastin Microfibril Interfacer 3	[52]
	FGFR3	Fibroblast growth factor receptor 3	[40]
	FOXC1	Forkhead Box C1	[81]
	GDF5	Growth differentiation factor 5	[82]
	ITGBL1	Integrin, beta-like 1 (with EGF-like)	[83]
	LECT1	Leukocyte cell derived chemotaxin 1	[84]
	TGFBR2	Transforming Growth Factor, Beta Receptor II	[85]
Muscle and tendon differentation	MYF5	Myogenic factor 5	[86]
	OLFM2	Neuronal Olfactomedin Related ER Localized Protein 2	[87]
	PLAC9	Placenta-specific 9	[88]
Axonal guidance and neural function	ACCN1	Acid-sensing ion channel 2	[89]
	CDH20	Cadherin 20, Type 2	[90]
	CPA6	Carboxypeptidase A6	[91]
	DAO	D-Amino-Acid Oxidase	[92]
	EPHA5	Ephrin type-A receptor 5	[93]
	ERMN	Ermin, ERM-Like Protein	[94]
	FAM19A5	Family With Sequence Similarity 19 (Chemokine (C-C Motif)-Like) Member 5	[95]
	FABP7	Fatty acid binding protein 7, brain	[96]
	GPR37L1	G Protein-Coupled Receptor 37 Like 1	[97]
	LPHN3	Latrophilin 3	[98]
	NRCAM	Neuronal cell adhesion molecule	[99]
	OTOS	Otospiralin	[100]
	PARK7	Parkinson protein 7	[101]
	RGMA	Repulsive Guidance Molecule Family Member A	[102]
	RGMB	RGM domain family, member B	[103]
	SLC7A10	Solute Carrier Family 7 (Neutral Amino Acid Transporter Light Chain, Asc System), Member 10	[104]
	TMEM47	Transmembrane Protein 47	[105]
Angiogenesis	ANGPTL1	Angiopoietin-related protein 1	[106]
	C1QTNF6	C1q and tumor necrosis factor related protein 6	[107]
	OSR1	Odd-skipped related 1	[108]
Regulation of tissue growth	CDKL1	Cyclin-dependent kinase-like 1	[109]
	CDKN2C	Cyclin-Dependent Kinase Inhibitor 2C	[110]
	LRIG3	Leucine-Rich Repeats And Immunoglobulin-Like Domains 3	[111]
	PDGFRL	Platelet-Derived Growth Factor Receptor-Like	[112]
	PTPRD	Protein tyrosine phosphatase, receptor type, D	[113]
	TMEM22	Solute carrier family 35 member G2	[114]
Genes with higher expression in EB and apical-LB tissues
Cell Signaling	ESRRG	Estrogen-related receptor gamma	[115]
	FAM150A	Family With Sequence Similarity 150, Member A	[116]
	FGF8	Fibroblast Growth Factor 8	[117]
	G3BP2	Ras GTPase-activating protein-binding protein 2	[22]
DNA Repair and Modification	HABP4	Intracellular hyaluronan-binding protein 4	[118]
	JARID2	Jumonji, AT Rich Interactive Domain 2	[62]
	OBFC2A	NABP1 nucleic acid binding protein 1	[119]
	ZMYM2	Zinc Finger, MYM-Type 2	[120]
Ubiquination	ASB6	Ankyrin repeat and SOCS box protein 6	[121]
	SENP1	SUMO1/sentrin specific peptidase 1	[122]
	SPSB4	SplA/Ryanodine Receptor Domain And SOCS Box Containing 4	[123]
	UBA3	Ubiquitin-like modifier activating enzyme 3	[124]
Transcriptional Regulation and RNA processing	HBP1	HMG-Box Transcription Factor 1	[125]
	NR2C1	Nuclear receptor subfamily 2, group C, member 1	[126]
	SFRS2	Serine/Arginine-Rich Splicing Factor 2	[127]
	SNRPA1	Small Nuclear Ribonucleoprotein Polypeptide A	[128]
	ZNF674	Zinc Finger Protein 674	[129]
Translation	ETF1	Eukaryotic Translation Termination Factor	[130]
	ISG20L2	Interferon stimulated exonuclease gene 20 kDa-like 2	[131]
	PUS3	Pseudouridylate Synthase 3	[132]
Mitochondrial Function	C21orf2	Chromosome 21 Open Reading Frame 2	[133]
	DNM1L	Dynamin-1-like protein	[134]
	TBC1D20	TBC1 Domain Family, Member 20	[135]
	WARS2	Tryptophanyl TRNA Synthetase 2, Mitochondrial	[136]
	YARS2	Tyrosyl-tRNA synthetase 2, mitochondrial	[137]
Nerve Related	MAPK6	Mitogen-Activated Protein Kinase 6	[138]
	MLLT11	Myeloid/Lymphoid Or Mixed-Lineage Leukemia (Trithorax Homolog, Drosophila); Translocated To, 11	[139]
	RTN4	Reticulon 4	[140]
	SMN1	Survival of motor neuron 1, telomeric	[141]

Discussion

ECM organization as a mechanism for positional plasticity and stability

We observed that the gene G3BP2 (Ras GTPase-activating protein-binding protein 2), which is part of a Twist1-G3BP2 mechanotransduction pathway, was expressed at higher levels in EB and apical-LB populations relative to basal-LB and stump populations. G3BP2 prevents Twist1 translocation to the nucleus, which results in activation of genes involved in differentiation [22]. Twist1 signaling is also regulated by matrix stiffness such that high stiffness in tissues results in G3BP2 release of Twist1, and activation of target genes. In the current study, we observed that the extracellular matrix molecule tenascin is more organized in the basal-LB tissue as compared to both EB and apical-LB tissue (Fig. 1); similar observations have been made for the blastema ECM as a whole [15]. It therefore is possible that the increased organization of the ECM in the basal-LB tissue alters Twist1-G3BP2 interactions, leading to an increase in Twist 1 nuclear translocation and expression of genes that promote differentiation in the blastema. Moreover, the increased abundance of genes involved in degrading the extracellular matrix (mmp1, mmp2, mmp3, mmp13, ctsk) that we have observed in EB and apical-LB could be upstream of G3BP2-mediated inhibition of differentiation in these plastic populations. Further experiments will be needed to test this hypothesis.

The intercalary response and positional plasticity

According to the Polar Coordinate Model of regeneration, the process of intercalation drives tissue growth. In our study we discovered that the amount of cell proliferation was higher in the EB and apical-LB populations relative to basal-LB and stump populations. Thus, we hypothesize that the intercalary response is occurring in the plastic blastema populations (i.e., EB and apical-LB). If this is true, then positional plasticity could be a consequence of intercalation. As we have explained previously, a number of genes that regulate epigenetic modifications and chromatin structure are enriched in plastic blastema populations. Given that positional information is epigenetically encoded [13, 23, 24], it is likely that the establishment of the position-specific epigenome in the newly generated cells is mediated by some of these genes and has a temporal component.

Our findings provide candidate genes and pathways that maybe involved in the regulation of cell proliferation. Our transcriptional analysis reveals that cndp2, which promotes cell proliferation, is more highly expressed in EB and apical-LB blastema populations. EGFR signaling is required for CNDP2 mediated proliferative activity [23] and is activated during regeneration [25]. Additionally, EGFR activity is induced by matrix metalloprotease activity in the regenerate [26]. We observed increased expression of a number of metalloproteases (mmp1, mmp2, mmp3, mmp13) in EB and apical-LB populations. These metalloproteases could be upstream of a CNDP2-mediated growth response during intercalation. Likewise, we also discovered that two growth inhibitory factors, cdkn2c and gas2 [27, 28], were more highly expressed in the basal-LB population. These results suggest the operation of mechanisms to inhibit growth in blastema tissues that are differentiating and not participating in a proliferation response.

Candidate genes and pathways for positional plasticity and stability

Several transcriptional studies of limb regeneration have been performed in recent years [29-33]. In all studies to date, mRNAs were isolated from partial or entire blastemas, including in some cases underlying mature stump tissue. While we did not achieve the temporal resolution of the Voss et al. 2015 analysis, our study is the first to compare global patterns of transcription among spatial domains of a blastema. Our results show that transcript abundances vary as a function of apical/basal position within the blastema. We also compared transcription among different blastema cell populations that were collected after performing proximal and distal amputation. Almost half of the genes were significant for the interaction term in our statistical model, thus suggesting that gene expression varied among blastema populations as a function of amputation site. If the proteins encoded by these genes contribute to the property of positional plasticity that is shared between EB and apical-LB cell populations, then this property may not correlate well with transcript abundance. However, we caution that many of the interaction gene effects were identified as significant because of our unbalanced experimental design – more replicates were collected and analyzed for proximal amputations, thus more statistical power was available to detect significant differences. Consistent with this interpretation, 324 post-hoc t-tests performed were significant when contrasting proximal cell population samples, while only 64 t-tests were significant when contrasting distal cell population samples. This pattern was also observed for genes identified as significant for the cell population term in the statistical model – 322 t-tests performed among proximal cell population samples were found to be significant compared to only 116 t-tests for distal cell population samples. Finally, only 24 genes were identified as significant for the amputation site term in the statistical model. While it seems likely that some genes maybe expressed differently between distal and proximal blastema cell populations, including genes that associate with positional plasticity, further studies with more replicate samples and finer proximal-distal and temporal sampling will be needed to rigorously test this idea. At any rate, it will be important in future transcriptional studies of limb regeneration to more finely investigate temporal and spatial components of variation.

While we consider all of 365 significant probesets identified from our statistical tests as potentially providing perspective about the molecular basis of positional information, we prioritized these genes further using t-tests. Working under the assumption that transcript abundance correlates with properties of positional information, we identified a high priority list of 65 genes. Highlighting genes primarily from this list, we discuss candidate genes and pathways that may underlie the properties of positional stability and plasticity.

Memory of positional information is a property that appears to be exclusive to fibroblasts in connective tissues throughout the body. This is based on transcriptional and epigenetic data [13, 23, 34], and the ability of connective tissues to induce the formation of new limb structures [10-12]. The EB-blastema mesenchyme is composed largely of cells of fibroblast origin [35] and thus the expression profiles in the EB mesenchyme should predominately reflect gene expression in blastema cells of fibroblast origin in a positionally plastic state. On the other hand, determining the transcriptional profiles that are linked with fibroblast-derived cells that have stabile positional information is more complicated because the basal-LB tissue is more heterogeneous, being composed of cells of multiple tissue origins [35-37]. However, since fibroblast-derived blastema cells contribute to the skeletal elements in the regenerate [36], the expression of genes involved in the development of cartilage/bone is likely associated with fibroblast-derived cells as they differentiate into these tissues. With this in mind, we discuss how genes that were significantly enriched in EB and apical-LB compared to basal-LB populations could control/affect the stability of positional information in fibroblast-derived blastema cells.

FGF signaling

FGF pathway members showed quantitative variation between positionally plastic and stable populations. At the EB stage, fgf8 expression initiates in basal cells of the wound epidermis and underlying blastema mesenchyme cells [38]. While fgf8 did not make the prioritized list of 65, expression was significantly higher in EB and apical-LB samples than basal-LB samples. This same pattern was observed for ebna1bp2, which encodes an FGF binding protein [39]. Inverse to these patterns, we observed higher expression of two FGF receptors (fgfr3, fgfrl1) and a binding protein that regulates FGF receptor signaling (fgfbp3) in the basal compartment of the LB blastema. These patterns suggest that FGF pathway components vary quantitatively between apical and basal compartments of a blastema. We hypothesize that these quantitative differences may contribute to a gradient mechanism that maintains positional plasticity in EB and apical-LB compartments, but promotes cellular differentiation and positional stability in the basal-LB compartment.

In support of this hypothesis, fgfr3 is expressed exclusively in differentiating chondrocytes in the proximal region of the developing mouse limb bud [40, 41]. One of the targets of FGF signaling during limb development is dact1, which we found to be highly expressed in basal-LB tissue. Dact1 is expressed in the limb bud mesenchyme, specifically in the developing cartilage [42-44]. Knockout of dact1 in mouse disrupts Planar Cell Polarity signaling [45], which plays a role in controlling the size and shape of the limb during limb bud morphogenesis [46]. FGFR3 and DACT1 colocalize to the same cells in developing limb bud, and FGF signaling activates the expression of dact1 [47]. SOX9, which is upstream of FGFR3 and COL9A2 expression and binds to cartilage-specific enhancers that regulate acan and col2A transcription [48], was also expressed more highly in basal-LB tissue. We observed an increase in the expression of both acan and col2A in the basal compartment, supporting the idea that similar mechanisms are underway in basal-LB blastema cells. It is possible that Sox9 activates the expression of fgfr3 and dact1 to induce the differentiation blastema cells into chondrocytes, and thus could also be linked to the stabilization of positional information in those cells during limb regeneration.

TGFb/BMP signaling

Given the role of TGFb/BMP signaling in the differentiation and development of the skeletal elements, it is not surprising to see an increased expression of a number of genes involved in this pathway in the basal compartment (tgfbr1, emilin3 crim1, gdf5 and osr1) where cartilage differentiation is actively underway. Tgfbr1, encodes for a TGFb receptor that plays an important role in skeletal development and regeneration in mammals and amphibians [49-51]. Emilin3 acts as a TGFb antagonist, and is expressed in the periskeletal tissue in developing limb buds [52]. Crim1 encodes for a gene that acts as a BMP antagonist by negatively regulating the processing and secretion of BMPs [53]. However, Crim1 is expressed in multiple tissues throughout the developing limb [54] and thus it is difficult to understand how this gene functions to regulate BMP signaling. Conversely, Gdf5 (aka BMP14) is a BMP-ligand that is expressed in, and is required for, developing and regenerating limb joints [55, 56]. Osr1 expression in the limb bud mesenchyme is essential for synovial joint formation, in part by maintaining the expression of gdf5 in the forming joint [57]. Interestingly, FGF signaling is required for osr1 expression in Xenopus embryos [58], and it is possible that some of the aforementioned FGF-related genes could play a role in the induction of osr1 expression in blastema cells as they differentiate into cartilage.

On the other hand, the inhibition of TGFb/BMP signaling could be required to prevent differentiation of positionally plastic cells. We observed higher expression of esrrg in the EB and apical-LB samples as compared to the basal-LB sample. This gene encodes an orphan nuclear receptor that binds to both the estrogen response element and the steroidogenic factor 1 response element, and activates gene expression in the absence of ligand [59]. ESRRG negatively regulates BMP-2 induced bone differentiation and formation in mammals [60], and thus this gene would be expected to function in preventing BMP-2 mediated differentiation of blastema cells.

Chromatin structure and epigenetic regulation

Large-scale chromatin rearrangements occur in cells of the uninjured limb as they dedifferentiate and become blastema cells, such that the chromatin becomes less densely packed and more euchromatic [61]. Consistent with this observation, genes that regulate epigenetic modifications (jarid2) and chromatin structure (zym2, habp4) are more highly expressed in EB and apical-LB tissues. Of note, JARID2 regulates differentiation of embryonic stem (ES) cells and cells within developing embryos by associating with PRC2 and inhibiting the trimethylation of histone H3 (H3K27Me3) [62]. JARID2 also plays a role in the recruitment of PRC2 to Hox loci in differentiating embryonic stem cells [63]. Presumably JARID2 plays a similar role in the blastema and is involved in the re-programming of positional information in the plastic blastema population. NFIA, which encodes for a transcription factor that can interact with and modify nucleosomal structure [64, 65], is enriched in the basal-LB tissue. This is interesting because NFIA promotes the differentiation of multiple tissue types [66-68], and likely has a similar function in the progenitor cells within the basal compartment. Altogether these observations indicate that genes that are involved with regulating chromatin structure and the epigenome are active in both plastic and stabilized blastema populations.

Ubiquination

Spsb4 and asb6 encode for proteins that are involved in protein ubiquination, and are more highly expressed in positionally plastic EB and apical-LB blastema populations. Since ubiquination is an important part of protein turnover [69], it is possible that the expression of these genes is linked to the increased metabolism of these dynamic cell populations (Fig. 2 and [70]). In addition, histone ubiquination is an important modification that alters nucleosome, and therefore chromatin structure [71]. Thus these ubiquitin–related genes may play a role in re-programming blastema cells as well as in cellular housekeeping.

Neural development and function

Above, we focused on genes that are likely to be expressed in fibroblast-derived blastema cells because they retain positional information. However, we note that a number of genes associated with neural development and function were expressed more highly in the basal-LB sample. These genes include asic2, cdh20, fabp7, fam19a5, llphn2, megf10, pcdh8, rgma, rgmb, scfd2, serpine2, tmem47. It is unclear whether these expression patterns reflect increased expression in the basal compartment, or a decreased expression the apical-LB compartment. However, if the abundance of these transcripts is indicative of neural activity in these tissues, these activities appear to be higher in basal-LB than apical-LB mesenchyme. We speculate that this could be a result of the formation of new neural connections with the differentiating tissues in this region.

Conclusion

In a previous study we discovered that positional information is unstable (or plastic) in EB and apical-LB populations, and is stabilized in basal-LB tissue [9]. To identify candidate mechanisms that regulate the stability of positional information in the blastema we compared EB and apical–LB samples to the basal-LB sample. We observed that genes associated with negative and positive regulation of cell differentiation were more abundantly expressed in plastic and stabilized populations, respectively. This suggests that the stability of positional information in blastema cells maybe associated with the regulation of signaling pathways for cellular differentiation. Our results also suggest that positional stability is a property of cells with increased ECM organization and associated gene expression. This observation likely reflects the increased amount of differentiating cells in the basal compartment, although it is also possible that ECM-linked mechanisms could regulate differentiation. From the current study it is unclear whether processes of cell differentiation and the stabilization of positional information are linked genetically, or whether the same molecular mechanism is utilized to “hardwire” both cellular and positional identity. We speculate that the latter possibility is likely, given the role of epigenetic modifications in the specification of both of these properties. Lastly, plastic blastema populations are associated with increased cell proliferation and higher expression of genes that promote growth. Since the theory of intercalation is based on the idea that positional discontinuity drives tissue growth, we posit that increased proliferation is an intercalary outcome. Altogether, our results provide an important first-step toward better characterizing spatiotemporal changes that occur in the developing blastema, and identifying candidate mechanisms that regulate the cellular property of positional-stability in the blastema during regeneration.

Methods

Animal husbandry

This study was carried out in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The experimental work was approved by the Institutional Animal Care and Use Committee of the University of California Irvine.

All of the experiments in this study were performed on the Mexican axolotl (Ambystoma mexicanum) measuring approximately 15–20 cm from snout to tail tip (Ambystoma Genetic Stock Center, UKY). Animals were anesthetized using a 0.1 % solution of MS222 (Ethyl 3-aminobenzoate methanesulfonate salt, Sigma), pH 7.0. To initiate regeneration, animals were either amputated just proximal to the carpals (distal amputation), or at the proximal end of the humerus (proximal amputation).

Tissue collection

The forelimbs and hindlimbs of animals were amputated at either the carpal level (distal), or at the upper arm (proximal). When the blastemas had reached early bud stage, the EB blastema mesenchymal tissue was harvested by gently teasing away the wound epithelium, and conservatively collecting blastema tissue, being careful not to collect mature limb tissue. On the same animals, the blastemas were allowed to regenerate until they reached LB stage. Similar to EB collection, the wound epithelium was carefully removed from LB blastemas. Blood vessels are sparse in the apical region of the LB, and the lack of blood vessels was used as an ocular cue for apical tissue collection. Basal-LB tissue was then collected conservatively as to not include stump tissue. See [9] for a more detailed description of how EB, apical-LB, and basal-LB tissue are collected. Stump tissue was collected from the same proximal/distal locations as the blastema tissues on unamputated animals, and included deep and superficial tissues in the limb. For each sample, tissues were pooled from 8 blastemas. Microarray was performed by the University of Kentucky microarray core using 3–4 independent biological replicate samples per each tissue type.

Microarray analysis

Samples were processed by RMA using Affymetrix Expression Console (Affymetrix, Santa Clara, CA). Samples were divided into three groups: early blastema (E) apical late blastema (A), and basal late blastema (B). All probesets that were consistently expressed below the bottom quartile of mean expression across all samples were removed from the analysis. Samples were also separated based on whether they were gathered from a distal (D) or proximal (P) amputation. Each probeset was fit to the generalized linear model of y = Βoc + Β1a + Β2 c * a + i + ε), where y was the measured log2 transformed intensity from the microarray, Β0, Β1, and Β2.were coefficients for the blastema cell type (c), amputation site (a), and blastema cell type by amputation site interaction effects, i was the intercept term, and ε was the normally distributed error. R was used to extract p-values from the full model using the lm function, as well as perform a two-way analysis of variance using the aov function to derive p-values associated with each main effect and the interaction term. Post-hoc t-tests to compare cell populations within an amputation site were also performed in R. The hclust and heatmap functions were used to cluster both samples and probesets using average linkage clustering on Pearson’s correlation [72]. Changes in gene expression were considered significant if the full model passed an FDR correction at 0.05, and post hoc t-tests were evaluated at p < 0.05. Probesets that were identified as showing a statistically significant change were clustered by average linkage clustering using 1-Pearson's correlation as a distance metric. Join distances were plotted as a function of the number of clusters to identify at what stage dissimilar expression profiles would be joined into the same group. From this plot, 15 clusters were selected and used [73, 74].

Gene enrichment analysis was performed using PANTHER gene expression tools [21] and the complete set of Biological Process gene ontology terms. The reference file for analysis was the complete list of annotated probesets (i.e., probesets with gene name acronyms) from the axolotl Affymetrix microarray. A familywise error rate (p < 0.00001) was used to correct for multiple tests.

Phalloidin and tenascin staining

Tissues for phalloidin-rhodamine and tenascin immunoflourescence were fixed, prepared for embedding in OCT, and cryosectioned into 10um slices. Sections were washed, and permeabilized using 0.01 % Triton-X prior to staining procedure. To stain the actin cytoskeleton, phalloidin-rhodamine (Cytoskeleton Inc., Denver, Co.) (14 uM) diluted 1:100 in PBS was incubated on sections for 1 h in the dark prior to washing and stabilizing with Vectashield mounting medium with DAPI (Vector Laboratories, Burlingham, Ca). Tenascin immunofluorescence was performed using the MT1 anti-newt-tenascin monoclonal antibody. The MT1 antibody, developed by Roy A. Tassava was obtained from the Developmental Studies Hybridoma Bank developed under the auspices of the NICHD and maintained by the University of Iowa, Department of Biology, Iowa City, IA 52242. Immunofluorescence was performed similar to [75] by incubating sections at 4 °C overnight with a 1:10 dilution in PBS of MT1 hybridoma media, washing, and a 1 h incubation with goat-anti-mouse antibody conjugated to Alexa-488 (Invitrogen). Images were obtained using a 20x objective on a Zeiss LSM780 (2-photon) confocal microscope.

Quantification of actin fiber characteristics

Confocal images were obtained of complete 10um deep, sagitally oriented, tissue sections of the blastema at different apical/basal positions, which were fluorescently stained for F-actin with phalloidin-rhodamine as described above. The tissue sections analyzed (Fig. 1b-c) were located approximately where the green lines indicate on Fig. 1a. The actin fiber orientations within each sagittal section were quantified using automated image processing implemented with a custom MATLAB script. The dominant orientation was estimated at each image pixel location by computing the image gradient using a Gaussian derivative filter (sigma = 2). The resulting orientation field was smoothed via standard methods [76] with a Gaussian filter (sigma = 10). To assess the differences in orientation across the spatial extent of each tissue section, the image was divided into small non-overlapping tiles (512x512 pixels) each covering approximately 150 square microns of tissue, and all together covered all of the mesenchymal tissue from each section. The distribution of gradient orientations in each tile was estimated using a histogram with 15 equal bins over the range of orientations (0–180°). Pixels with a gradient magnitude of less that 1 % of the maximum image intensity were excluded from analysis. This threshold served to eliminate locations in the image where no stained actin fibers were visible by eye (and where the corresponding orientation estimate is unstable). The discrete entropy of each histogram was computed as a summary statistic to measure the degree of order (alignment) or disorder of actin filaments within the region of tissue spanned by the tile.

The average of the discrete entropy was calculated for each apical/basal region of the blastema by calculating the average entropy of all of the tiles that were of only mesenchymal tissue (i.e., did not include epithelium) from a tissue section at each location. Figure 1c is a representative data set of the quantification of the entropy values from individual tissue sections at different positions within a single blastema. The quantification of the actin cytoskeleton was performed on multiple blastemas that were sectioned either sagitally or transversally, and each showed the same trend of decreasing entropy in the increasingly basal tissues. The data from multiple blastemas was not combined because we could not ensure that the apical/basal location of each section was equivalent from blastema to blastema. In Additional file 5: Figure S2 we provide an example of the workflow on the analysis that was performed on a transverse tissue section of a late blastema. Note that the same apical/basal trend of actin organization was observed.

Edu injection, detection, and quantification of labeling index

100 ug of Edu (Invitrogen, Carlsbad, Ca) was injected intraperitoneal in the flank of the animal 3 h prior to tissue collection. Tissues were fixed in 3.7 % PFA solution, embedded in paraffin, and detected using the manufactures protocol. Vectashield with DAPI (Vector Laboratories, Burlingham, Ca) was used to stain all nuclei. A tile scan of the entire blastema was obtained using a 20x objective on a Zeiss LSM780 (2-photon) confocal microscope. To quantify the % labeling index of the blastema mesenchyme, the wound epithelium was removed digitally using Photoshop, and the blastema was divided into digital sections approximately 200 uM along the proximal/distal axis. Edu and DAPI stained nuclei were manually counted using the cell counter in the ‘analyze’ plug-in on ImageJ. The % labeling index was determined for each digital section by calculating the ratio of Edu positive cells to total number of cells, and multiplying by 100.

Availability of supporting data

The microarray data files supporting the results of this article are available for download at Sal-Site http://www.ambystoma.org/genome-resources/20-gene-expression.