Retinoic acid receptor regulation of epimorphic and homeostatic regeneration in the axolotl

Urodele amphibians (salamanders) are capable of regenerating amputated limbs and tails throughout life by recruiting cells juxtaposed to the amputation plane to migrate distally (towards the hand) and proliferate into a mass of lineage-restricted cells called a blastema (Kragl et al., 2009; Monaghan and Maden, 2012a). Blastemas only regenerate structures distal to their origin, known as the ‘rule of distal transformation’, using positional cues provided by cells proximal to the amputation plane (Ludolph et al., 1990; Maden, 1980; Stocum and Thoms, 1984). Young blastemal cells are in a state of cellular plasticity, which allows them to adopt distal positional values (McCusker et al., 2014; McCusker and Gardiner, 2013; Roensch et al., 2013). Young distal limb blastema cells can be reprogrammed with supplemental retinoic acid (RA) to a proximal fate (Maden, 1982), posterior fate (Kim and Stocum, 1986; Stocum and Thoms, 1984) and ventral fate (Ludolph et al., 1990), which will not occur in uninjured limbs (McCusker et al., 2014; Niazi et al., 1985) or after redifferentiation has commenced (Niazi et al., 1985). Despite the power of this experimental approach for understanding the role of RA during regeneration and how positional identity is established and maintained, little is known about the transcriptional network that regulates positional information.

RA is a molecule with pleiotropic functions that is vital during vertebrate development for regulating embryo patterning, cell differentiation, and organogenesis (Clagett-Dame and DeLuca, 2002; Duester, 2013; Marlétaz et al., 2006). RA signaling controls developmental processes by regulating gene transcription through the activation of retinoic acid receptors (RARα, RARβ and RARγ). RARs heterodimerize to retinoid X receptors (RXRs) and, together, these transcriptional complexes bind to retinoic acid DNA response elements (RAREs) located near RA target genes (Chambon, 1996). RAR/RXR complexes work as transcriptional repressors with no ligand and as activators in the presence of RA ligand (Rochette-Egly and Germain, 2009). Limiting RA concentration, inhibiting RAR signaling, or inhibiting RA metabolism has detrimental effects on limb development in chicks (Roselló-Díez et al., 2011; Stratford et al., 1996), zebrafish (Grandel et al., 2002) and mammals (Dranse et al., 2011; Lohnes et al., 1994; Niederreither et al., 2002; Sandell et al., 2012, 2007; Williams et al., 2009; Yashiro et al., 2004). The role of RA during limb regeneration is less clear (Blum and Begemann, 2013), although several lines of evidence support an active role. RA is present in regenerating limbs (Scadding and Maden, 1994) and RA-reporter axolotls show RA signaling in regenerating limbs (Monaghan and Maden, 2012b). Genes that regulate RA signaling are expressed in regenerating frog limbs (McEwan et al., 2011) and salamanders including Rdh10 (Monaghan et al., 2012), Raldh1 (Knapp et al., 2013), Raldh3 (Monaghan et al., 2012), Rarα (Ragsdale et al., 1989) Rarβ (Carter et al., 2011; Giguère et al., 1989) and Rarγ (Hill et al., 1993; Ragsdale et al., 1989; Voss et al., 2015). Also, Raldh inhibition blocks axolotl limb regeneration (Maden, 1998; Scadding, 2000), and also epimorphic fin zebrafish regeneration (Blum and Begemann, 2012), and excess RA induces duplicated patterning during Xenopus hindlimb regeneration (Cuervo and Chimal-Monroy, 2013) as it does in salamanders.

RA will reprogram regenerating limbs up to the early/mid limb bud stage in axolotl salamanders, but generates hypomorphic limbs when treated during development. Both phenotypes can be observed simultaneously in axolotls because hindlimbs emerge late in development, when forelimbs have already completely differentiated (Scadding and Maden, 1986). Hypomorphic regeneration also occurs when RA is added to limbs after the early/mid bud stage, suggesting that RA signaling cannot influence the PD axis outside the developmentally plastic phase of the early blastema (Maden, 1983; Niazi et al., 1985). Our previous work using reporter-based analysis supports this hypothesis because RA reporter activity is different between developing and regenerating limbs. Furthermore, adding excess RA during the early bud stage of regeneration (5 days post amputation) induced RA reporter activity in blastema connective tissue fibroblasts (Monaghan and Maden, 2012b), the precise cells responsible for PD duplications (Nacu et al., 2013). Similar to the effects of adding RA after the early/mid bud stage has commenced, Rarβ antagonism with the isoform-specific antagonist LE135 has no effect in early regeneration, but halts regeneration at the mid/late bud stage (Del Rincón and Scadding, 2002). Therefore, the differential effect of RA on developing and regenerating limbs might be due to its interactions with specific RARs during specific stages of regeneration or in specific cell types. RA's teratogenic capacity to truncate limbs rather than re-specify PD axis identity could be explained by dysregulation of specific RARs. It is fundamental to our understanding of limb development and regeneration to identify the molecular basis of proximodistal duplication versus truncation of the regenerating limb.

The cellular mechanisms that impart positional memory are still unclear (McCusker et al., 2015; Phan et al., 2015; Roensch et al., 2013). Several transcription factors have been identified that presumably activate genes responsible for positional memory (Crawford and Stocum, 1988), including Meis1, Meis2 (Mercader et al., 2005), Hoxd10 (Simon and Tabin, 1993) and Hoxa13 (Gardiner et al., 1995), but our understanding of what makes a limb proximodistally duplicate, truncate, or grow the proper structure is lacking. Fundamental questions are unresolved including how many genes participate in PD positional memory, how these genes are coordinated at the cellular level, and whether salamander orphan genes regulate the positional memory required for regeneration. Thus, the objective of this study was to reveal the underlying basis of RA-induced PD duplications versus truncations utilizing transcriptomics, RARE-reporter animals, and RAR-specific agonists and antagonists.

Our previous work showed that RAR reporter activity is present in regenerating limbs with expression mainly in epidermal keratinocytes, axons and nerve-associated cells. RA-induced PD duplication coincided with upregulation of RA reporter activity in connective tissue fibroblasts (Monaghan and Maden, 2012b). Here, we investigated whether endogenous RAR activity is required for limb regeneration. We treated regenerating animals with the selective RARβ antagonist LE135 (Li et al., 1999), because it has been shown to cause limb truncations and hypomorphic regenerates whereas RARα-specific and pan-RAR antagonists have minimal effects (Del Rincón and Scadding, 2002). At 7 days post amputation (dpa), RARβ antagonism induced reporter activity in RARE-EGFP limbs to a similar extent as RA-treated limbs, rather than decreasing activity as would be expected [Fig. 1A-C; n=6, ∼4 cm total length (TL)]. Reporter activity was mainly present within skeletal tissue including the perichondrium, in a few fibroblast-like cells, and within the basal wound epidermis compared with basal wound epidermis and fibroblasts in RA-treated limbs (Fig. 1C). Overall, RARE reporter activity had similar patterns of expression as RA-treated limbs except that LE135 induced RARE-EGFP more strongly in skeletal tissue.

RARβ antagonism did not halt blastema formation or initial growth. Rather, LE135 significantly halted growth at the mid bud blastema stage (Fig. 1D-F′), which corresponds approximately to the beginning of re-differentiation. After 15 days of treatment with a different RARβ antagonist, LE540 (Li et al., 1999), blastema size was 1.044±0.16 s.d. (n=5 right limbs) versus 1.312±0.12 s.d. in untreated limbs (n=6 right limbs) (Student's t-test, two-tailed; P=0.01) and had progressed to pallet stage compared with early digit formation in untreated limbs. Based upon these observations, we reasoned that RARβ inhibition might negatively impact endochondral ossification. Alcian Blue staining showed that some cartilaginous precursors (chondroblasts) are formed in LE135-treated limbs (Fig. 1E′ versus1F′) along with the expression of Collagen 2a protein (Fig. 1G versus 1H), but LE135-treated limbs showed a lack of progression from chondroblasts to chondrocytes as indicated by the formation of lacunae-like structures as shown in Fig. 1G (arrowheads) versus Fig. 1H.

We next tested whether RARβ antagonism also inhibits limb development through an RA-responsive transcriptional pathway. We found that RARβ antagonism, initiated at the onset of forelimb bud outgrowth (stage 36), slowed forelimb growth by the mid bud stage 43 (Fig. 1I) and growth ceased by stage 50 (Fig. 1J,K). RA reporter activity is known to be present in developing forelimbs, but is absent in developing hindlimbs (Monaghan and Maden, 2012b). We found that RARβ antagonism initiated at the onset of hindlimb bud outgrowth (stage 51) activated RARE-reporter activity 6 days later (Fig. 1L), and this corresponded with inhibition of hindlimb outgrowth (Fig. 1M). Reporter activity was increased in the epidermis and proximal mesenchyme, which is the region of chondrocyte differentiation in the developing limb (Fig. 1L). This shows that despite RA signaling having different in vivo patterns between forelimb and hindlimb development, RARβ antagonism generates a similar outcome. Several mechanisms might explain the induction of RARE activity after RARβ antagonism. One possibility is that LE135 is acting as an RARβ agonist instead of antagonist. A second possibility is that RARβ has an inhibitory role in the absence of ligand, normally preventing transcription of target genes, but when this inhibitory activity is inactivated, gene expression of RA target genes is induced. A similar inhibitory role of RARs occurs during mammalian chondrogenesis, when adding RAR antagonists induces gene expression of some RAR target genes (Weston et al., 2002, 2003b). Therefore, it is possible that transcriptional inhibition was removed with RARβ antagonism, inducing RARE-dependent gene expression programs.

To test whether RARβ antagonism induces RARE-dependent transcriptional changes as well as delineate the molecular basis of RA-induced PD duplication versus truncation, we performed microarray gene expression analysis on forelimbs that will eventually regenerate normally (DMSO treated), become PD duplicated (RA treated) or become truncated (LE135 treated; Fig. 2A). Genes were identified as statistically significant if they had a false discovery rate (FDR)<0.05 determined by an ANOVA analysis (533 significant probe sets), and exhibited a >1.5-fold change (FC) relative to control DMSO samples in either treatment group (327 significantly changed probe sets). Surprisingly, high similarity was observed in gene expression between LE135- and RA-treated forelimbs despite yielding different phenotypes (Pearson's correlation coefficient between treatment groups using log2 fold change from DMSO=0.883). Pairwise comparisons between groups (FC>1.5 and FDR<0.05) showed that most genes upregulated after RA treatment were also upregulated after LE135 treatment (Fig. 2B) suggesting a similar transcriptional ‘activating’ response in both treatment groups. Many more genes were uniquely downregulated between RA- and LE135-treated groups suggesting a more divergent transcriptional ‘silencing’ response between PD duplication and truncation (Fig. 2C).

To classify quantitative differences between treatments, hierarchical clustering of significant genes was performed on all 327 significantly changed genes, which generated five distinct clusters (Fig. 2D). Cluster 1 (n=97) were on average upregulated after both treatments compared with controls. Cluster 2 (n=67) included genes that were on average higher in RA-treated samples. Cluster 3 (n=34) included genes that on average were unchanged in LE135-treated samples, but were upregulated after RA treatment (Table S1). In contrast, cluster 4 (n=14) contained genes that on average changed little after RARβ antagonism, but were downregulated during RA-induced PD duplication (Table S1). Lastly, cluster 5 (n=115) contained genes that were on average downregulated in both treatment groups. Overall, hierarchical clustering highlighted the dynamic transcriptional response that occurs after perturbation of RAR signaling.

We will first focus on common gene expression changes observed after either treatment. The most strikingly upregulated genes in both treatment groups were involved in the retinoic metabolic process (over-representation analysis) including genes involved in RA synthesis, shuttling to the nucleus, catabolism, and RA-dependent transcriptional activation and repression (Fig. 2E). This suggests that RA signaling increases in both treatment groups, even though RA was not introduced to LE135-treated limbs. One explanation for this is the upregulation of RA synthesis genes after LE135 treatment (Fig. 2E). Another group of commonly upregulated genes were involved in sterol metabolism including Cyb5a, Soat1, Sdr16c5, Dhrs3, Gmds and Cyp4b1. Other striking expression patterns in cluster 1 included the upregulation of genes associated with extracellular matrix production and breakdown including Aggrecan, Brevican, Efemp1, Elfn1 and Mmp13 as well as intracellular intermediate filaments including Krt8, Krt15 and Krt19. It is clear that some common cellular changes are occurring in both PD duplicated and truncated limbs.

We reasoned that identifying genes specifically induced or silenced during PD duplication compared with controls would reveal the underlying mechanism of RA-induced PD duplication. Therefore, we focused on clusters 2-4, which included differentially regulated genes in RA-treated limbs compared with LE135-treated and DMSO-treated limbs. Although LE135 may have upregulated some of the same genes, clusters 2 and 3 show that the level of upregulation is on average much lower than RA-treated limbs. This might be due to the fact that many RA synthesis genes are upregulated after LE135 treatment. Many cluster 2/3 genes (n=101) are expressed in proximal developing limb buds in other limbed vertebrates or required for proper limb development [cluster 2 expressed in proximal limb: Meis1 (Mercader et al., 2000, 2005), Meis2 (Mercader et al., 2000, 2005), Pbx1 (Selleri et al., 2001), Arid5b (Ristevski et al., 2001); cluster 2 expressed in limb bud: Mia3 (Bosserhoff et al., 2004), Rac1 (Bell et al., 2004; Suzuki et al., 2013), Asph (Patel et al., 2014), Neo1 (Hong et al., 2012), Cyp26B1 (MacLean et al., 2001), Flrt2 (Haines et al., 2006), Rarγ (Pennimpede et al., 2010), Rbp1 (Gustafson et al., 1993), KIAA1217 (Semba et al., 2006); cluster 3 (Table S1) expressed in proximal limb: Fibin (Taher et al., 2011; Wakahara et al., 2007), Epha7 (Araujo et al., 1998), Nrip1 (Smith et al., 2014), Rnd3 (Bell et al., 2004); cluster 3 expressed in limb bud: Apcdd1 (Jukkola et al., 2004), Zfn638 (Bell et al., 2004), Stat3 (Gray et al., 2004), Tsh2 (Caubit et al., 2000; Erkner et al., 1999)]. The association of these genes with limb patterning in other vertebrates supports the idea that RA reprograms the distal cells to resemble a proximal limb cell fate. It also suggests that PD duplication entails at least 100 genes. Genes that have been previously identified as upregulated after RA treatment in regenerating salamander limbs were also identified in our study including Meis1 and Meis2 (Mercader et al., 2005; Simon and Tabin, 1993), genes that are accepted as determining proximal fates in vertebrate limbs (Mercader et al., 2000; Roselló-Díez et al., 2011) (Meis1 FC=+1.99 RA, FC=+1.35 LE135; Meis2 FC=+1.92 RA, FC=+1.52 LE135).

Cluster 4 included 14 downregulated genes in RA-treated limbs compared with LE135-treated and DMSO-treated limbs (Table S1). Alox5 was the only exception because it was exclusively upregulated in RARβ antagonized limbs (LE135 versus DMSO +1.55-fold; RA versus DMSO −1.17). Seven of the 13 genes downregulated in RA-treated limbs are known to be expressed in the distal portion of the developing or regenerating vertebrate limb including Lhx9 (Gu and Kania, 2010; Tzchori et al., 2009), Zic5 (Merzdorf, 2007), Lmo1 (Taher et al., 2011), Lhx2 (Taher et al., 2011; Tzchori et al., 2009), Spry1 (Minowada et al., 1999; Wang and Beck, 2014), Msx2 (Bell et al., 2003; Carlson et al., 1998; Tribioli et al., 2002), HoxA13 (Gardiner et al., 1995; Haack and Gruss, 1993; Scotti et al., 2015), most of which are required for distal identity in developing mouse limbs. This suggests that distal-identity genes are silenced only in limbs undergoing PD duplication, similar to the transcriptional activation of proximal-identity genes during PD duplication.

Positional information is thought to reside on the cell surface of blastema cells (Stocum and Cameron, 2011) or in the extracellular matrix (Phan et al., 2015), which is supported by the fact that proximal blastemas engulf distal blastemas in vitro (Nardi and Stocum, 1984). Our data provide several candidate molecules for regulating positional information in clusters 2, 3 and 4 (n=115), which included 23 extracellular molecules (GO term: Extracellular Region) as well as 11 genes involved in the regulation of cell adhesion (GO term: Cell Adhesion). Overall, microarray analysis supports the idea that PD duplication entails both loss of distal cell identity and gain of proximal cell identity, and modifications in cell-cell contact and cell adhesion properties.

Rarγ in particular has been associated with regulating PD limb duplications (Pecorino et al., 1996). Our results show that RA-induced PD duplication increased Rarγ expression to 1.62-fold higher than controls (cluster 3) versus 1.30-fold in LE135-treated limbs. qPCR supports this finding and shows that RARα and RARβ are not upregulated in either treatment group (Fig. 2G). Previous work has shown that activation of RARδ alone, which is homologous to human RARγ, was able to proximalize cells whereas RARα and RARβ were incapable (Pecorino et al., 1996). To test whether activation of RARγ is also capable of proximalizing entire limb blastemas, we treated early blastemas with a potent RARγ selective agonist, CD1530. We find that RARγ agonist treatment of early limb blastemas was capable of mimicking RA treatment by generating PD duplications to the shoulder level (n=2; Fig. 2H). This result supports the hypothesis that RARγ is the key RAR regulating the PD limb axis during limb regeneration, although a more thorough analysis of other RAR agonists and antagonists is clearly needed to support this claim.

RARβ antagonism inhibited limb growth leading to limb truncation during development and regeneration. Genes associated with limb truncation were found mainly in cluster 1 (Table S1). The first striking feature of cluster 1 is that it contains genes involved in skeletal formation and remodeling including the osteoblast master regulator gene Sp7, which is higher after RARβ antagonism (FC=+1.77 after LE135 treatment versus FC=+1.07 after RA treatment). Other genes known to be upregulated after osteoclastogenesis included tank (Maruyama et al., 2012) (FC=+1.51 after LE135 treatment), and lipid mediators including Alox5 (cluster 5), Alox15b and Aloxe3. In mammals, loss of lipid mediators Alox5 and Alox15b leads to an increase in bone, and increase of the activity of these lipid mediators decreases bone density (O'Connor et al., 2014). Other gene expression patterns were suggestive for an effect on skeletal progenitor differentiation including a FC of +2.2 of Tgfβ2 in LE135-treated limbs (FC=+1.73 in RA treated), a FC of +1.47 of Tgfβ1 in LE135-treated limbs (FC=+3.30 in RA treated), and a significant downregulation of Bmpr1b (FC=−2.57 in LE135 and FC=−1.77 in RA-treated limbs). Although most differences between RARβ antagonism and RA-induced PD duplication were quantitative in nature, it seems that gene expression patterns were skewed towards a transcriptional program leading to skeletal regression.

Considering the lack of skeletal differentiation that occurs in regenerating limbs after RARβ antagonism, we next investigated whether RARβ antagonism has an impact on uninjured bone integrity. RARβ antagonism led to a permanent shrunken limb phenotype (Fig. 3A-C). After 21 days of treatment in smaller animals, severe shrinking occurred [n=8 controls, snout to vent length (SVL)=2.5, TL=4.8, control stylopod+autopod=5.33±0.58 s.d., treated stylopod+autopod=2.49±0.82 s.d.; Student's t-test, two-tailed; P<0.001]. Integrity of long bones was strikingly impacted compared with untreated limbs, which was associated with an increase in RARE-EGFP reporter activity in long bones, epidermis and nerve axons (Fig. 3D), suggesting that the effect of LE135 could be partially cell intrinsic. Effects of RARβ antagonism included a compaction of the metaphysis and diaphysis with little effect on the epiphysis and an increase in osteoclasts within the diaphysis of bones (Fig. 3F,G). Defects were clearly apparent after microCT evaluation at 12 days of treatment, although no significant decrease in radius/ulna length could be observed at this point. Overall, the loss of bone homeostasis is consistent with gene expression profiles described in the results section above. Animals had an excess amount of skin, suggesting that degeneration was specific to the skeleton (Fig. 3H versus 3I). Furthermore, the cartilaginous epiphysis of treated limbs and carpals of the hands were of normal size (Fig. 3A versus 3B) suggesting degeneration of differentiated chondrocytes. Overall, long bone degeneration caused by RARβ antagonism seems to be due to an active transcriptional response within the differentiating skeletal cells, which is associated with significant osteoclastogenesis.

We next investigated whether the negative impact of RAR perturbation was specific to the limb. LE135 treatment for 21 days resulted in scoliosis of the spine demonstrating that the effects of RARβ antagonism also occurred in other skeletal tissues (Fig. 4A-C). RARE-EGFP animals show that RA reporter activity is minimal in the uninjured spinal column, except in spinal cord axons and a few cartilage cells (Fig. 4D). Upon RARβ antagonism, reporter activity increased in chondrocytes surrounding the spinal cord, especially in the dorsal chondrification center of the neural arch (Fig. 4E). In contrast, RA treatment induced reporter activity primarily in neural progenitor cells of the spinal cord, some white matter cells, the neural meninges, and cells resembling fibroblasts in the muscle (Fig. 4F). Altogether, these data strongly suggest that RARβ antagonism induces a specific RA-transcriptional response in skeletal tissue, which leads to a loss of skeletal integrity, possibly through a loss of homeostatic regenerative ability. RA induces a more specific response in fibroblastic cells, supporting the idea that RA specifically reprograms fibroblast cell identity.

Based upon the similar RAR-dependent reporter activity in uninjured tails and limbs, we next assessed whether RARβ antagonism also impacts tail regeneration. Indeed, RA reporter activity was primarily localized in axons of the untreated regenerating spinal cord (Fig. 4G-I), whereas RARβ antagonism induced significant reporter activity in differentiating prechondrocytes and epidermis (Fig. 4J-L). RA treatment increased reporter activity in spinal cord neural progenitor cells and some fibroblasts (Fig. 4M-O), which could explain the inhibitory properties of RA on spinal cord cell proliferation and urodele tail regeneration (Pietsch, 1993). RA is also known to regulate neural differentiation across vertebrates (Maden, 2007). The similar responses of RA treatment and RARβ antagonism between the limb and tail suggests that there may be a common RAR gene expression program regulating both limb and tail regeneration. Overall, the contrasting cell types responding to RA treatment versus RARβ antagonism also suggests that the role of RARs during regeneration is partially cell type dependent.

In this study, we show that modulation of RAR activity has a significant impact on tissue patterning and differentiation during epimorphic regeneration and skeletal homeostasis. We utilized reporter animals and gene microarrays to show that pharmacological activation of RARs with RA treatment, presumably through RARγ activation (Fig. 2G), induced a proximalization program leading to limb PD duplications. RARβ antagonism negatively affected skeletal differentiation and growth during epimorphic limb and tail regeneration and induced a skeletal regression program in uninjured skeleton. RARE-EGFP animals showed that induction of each transcriptional program had some overlap between tissue types, but also showed unique expression patterns – chondrocytes in the case of the truncation program and fibroblasts in the case of the PD duplication program (Fig. 4P). We propose that proper RAR activation is essential in a cell type-dependent and temporal manner. Overall, highly regulated RAR activity controls crucial transcriptional networks required for tissue patterning, differentiation, and tissue turnover during both epimorphic regeneration and homeostasis.

The endogenous role of RARs during tissue regeneration is unclear. We show that an RARγ agonist alone is sufficient for PD limb duplications, suggesting that RARγ might regulate patterning. This is supported by a microarray study (Voss et al., 2015) showing that Rarγ transcripts increase at the onset of blastema formation and stabilize thereafter. qPCR analysis also shows that only RARγ, not RARα or RARβ, is upregulated during PD duplication (Fig. 2G). These results together reinforce findings that RARγ is capable of proximalizing distal newt blastema cells, but RARα and RARβ cannot (Pecorino et al., 1996), and the fact that RARα antagonists have little impact on axolotl limb regeneration (Del Rincón and Scadding, 2002). It is possible that RARγ activity sets the appropriate PD level of the early blastema and overactivation with agonists sets the level to a proximal fate.

Few studies have screened for genes involved in positional re-specification of the limb. One exception used subtractive cDNA screening to identify upregulated and downregulated genes in distal newt blastemas after RA treatment (da Silva et al., 2002). This study identified one salamander-specific molecule (Geng et al., 2015), prod1, that has a PD gradient in newts and can proximalize distal blastemal cells in newts and axolotls (Echeverri and Tanaka, 2005). We did not observe an upregulation of axolotl Prod1 after RA treatment, which supports the finding that Prod1 transcripts are more abundant in distal blastemas compared with proximal blastemas in axolotls (McCusker et al., 2015). The current model is that Prod1 signals through epidermal growth factor receptor to induce Mmp9 expression (Blassberg et al., 2011). In our study, Mmp9 was not differentially regulated between treatment groups although it is upregulated during the early stages of limb regeneration (Monaghan et al., 2009; Yang et al., 1999). Considering that Prod1 is predicted to be a secreted molecule in all other salamanders (Blassberg et al., 2011), it will be important to test whether it plays an endogenous functional role in the axolotl and is required for PD limb patterning as it is in newts (Kumar et al., 2015).

One model for vertebrate limb patterning is that trunk-derived mesoderm generates a proximal source of RA, which induces expression of the stylopod-specific homeobox genes Meis1 and Meis2 (Cooper et al., 2011; Rosello-Diez et al., 2014; Roselló-Díez et al., 2011). RA signaling is inhibited distally by Fgfs (Cooper et al., 2011; Mariani et al., 2008) and Cyp26b (Yashiro et al., 2004), which is supported by genetic ablation of distal Fgf genes (Mariani et al., 2008) or Cyp26b (Yashiro et al., 2004). Our data partially support this model as we observed clear upregulation of proximal Meis1 and Meis2 genes and the downregulation of Sprouty1, a gene upregulated by FGF signaling after RA treatment (Minowada et al., 1999; Wang and Beck, 2014). Furthermore, clusters 2-4 clearly showed an induction of proximally expressed genes and silencing of distally expressed genes. The permanent change in PD cell identity is likely to require restructuring of the epigenetic landscape. In support of this hypothesis, we found that Ncoa3, the key ligand-dependent co-activator of RAR target genes (Torchia et al., 1997), was upregulated during PD duplication (Fig. 2E,F). Furthermore, Nrip1, the key ligand-dependent co-repressor of RAR target genes (Hu et al., 2004) was also upregulated (Fig. 2E,F) as well as the downregulation of the histone methyltransferase Whsc1 (Nimura et al., 2009) and differential expression of many homeobox-containing genes (Meis1, Meis2, Pbx1, Hoxc5, Tshz2, Zhx1, Zfhx4, Msx2, Hoxa13, Lhx2, Dlx6 and Lhx9). Together, this group of genes is likely to be crucial for re-specification of positional information in the regenerating limb.

The similarity in gene expression between PD duplication and truncation (Fig. 2) was surprising considering the divergent phenotypes. For example, genes associated with the proximal identity of vertebrate limbs, including Meis1, Meis2 and Pbx1, were upregulated in both treatment groups. This may be explained by the fact that RA synthesis genes are upregulated after LE135 treatment and Meis expression is due to new RA synthesis. Alternatively, it could be accounted for by the fact that Meis proteins are expressed after axolotl limb amputation in muscle blastema cells (Nacu et al., 2013) and epithelium (Nacu et al., 2016), which probably respond differently than fibroblast-expressing Meis. Another possible scenario is that RARβ antagonism might partially reprogram PD identity, but the program is incomplete or the truncation transcriptional program overrides the PD program. Regardless, genes found in cluster 3 including Tshz2, Tll2, Htra3, Fibin and Cetp might be new indicators of limb proximalization, supplementing classical indicators of proximal limb identity. A limitation of our study is that whole blastemas were analyzed rather than fibroblasts specifically. It would be interesting in the future to assess global gene expression changes only in fibroblasts, which are the cells known to regulate positional information of the limb.

Our data suggest that the mechanism by which LE135 inhibits epimorphic regeneration is through disruption of endochondral ossification. This leads to the question of how an antagonist can increase RAR target gene expression. During chondrogenesis, RARs play a repressive function; ligand-less RARs/RXRs recruit repressive transcriptional complexes to RA target gene promoters, which allow the chondrogenesis program to progress. In vitro, RAR-mediated repression is required for chondrocyte differentiation (Weston et al., 2003a, 2002). Chondrogenesis is also inhibited by agonists for RARα (Shimono et al., 2010; Weston et al., 2002) or RARγ (Shimono et al., 2011; Williams et al., 2009) (promotes RAR transcriptional activity) and enhanced by RAR reverse agonists (Williams et al., 2009) (promotes RAR transcriptional repression). In Cyp26b1 null mice (excess RA), skeletal prechondrocytes begin to differentiate, but exhibit reduced chondrocyte differentiation (Dranse et al., 2011). In our study, a similar mechanism might occur in that LE135 inhibits the repressive function of RARβ, activating the wrong transcriptional program in prechondrocytes (cluster 1 and Alox5). This could account for the similar gene expression patterns observed between RA treatment and LE135 treatment.

In vertebrates, long bones undergo continuous turnover, otherwise known as homeostatic regeneration, through osteoblast-based addition and osteoclast resorption. Excessive RA signaling is known to impact homeostatic turnover and skeletal integrity of long bones, including conditions like hypervitaminosis A (Green et al., 2016; Henning et al., 2015). Excess RA signaling increases osteoclast formation in mammals in vitro and in vivo (Henning et al., 2015); this is also observed in our studies after LE135 treatment i.e. increased RA reporter activity in skeletal tissue (Fig. 3D and Fig. 4P), increased osteoclastic gene expression, and increased numbers of osteoclasts in resorbing bone (Fig. 3F,G). Furthermore, in vivo data suggest that loss of RAR repression leads to accelerated chondrocyte hypertrophy (Dranse et al., 2011), which we also observed after LE135 treatment (Fig. 3F). It is likely that in our studies, increased RA signaling is context dependent – RA ligand-based RA signaling might not shrink skeletal tissue, but LE135-induced transcription does promote resorption. The transcriptional responses specific to LE135 treatments should provide insight into RA signaling-induced bone resorption.

Our study further elucidates the roles of RARs during regeneration, but also brings to light several unknowns about limb regeneration. The most pressing of which is whether endogenous RA ligands are required for limb regeneration and whether the PD duplication of the limb is exclusively regulated by RARγ. Furthermore, it will be important to determine the functions of genes regulated by RARs during PD duplication; are they capable of determining proximodistal identity and are they required for the process? The results presented here provide crucial information for tackling these problems.

Ambystoma mexicanum (axolotls) were bred in captivity either at the University of Florida or Northeastern University. Experiments were performed in accordance with University of Florida and Northeastern University Institutional Animal Care and Use Committees. For all experiments, animals were anesthetized by treatment of 0.01% benzocaine. In all cases of amputations, the radius/ulna or femur were trimmed to make a flush amputation plane and limb staging was performed according to Armstrong and Malacinski (1989) and Nye et al. (2003). Animals were bathed in drug [RA, 1 μM (Sigma); LE135 (Tocris), 250 nM; CD1530 (Tocris), 250 nM; LE540 (Wako), 1 μM; 0.03% DMSO (Sigma)] for the designated times with water changes every other day or every day for the microarray experiment.

RARE-EGFP sections were fixed in 4% paraformaldehyde at 4°C overnight, cryomounted in OCT medium (TissueTek), sectioned at 15-20 µm, stained in Hoechst 33258, and mounted in 80% glycerol. Histology was performed by fixing tissues in 10% neutral buffered formalin at 4°C overnight, washing twice in PBS, processing for paraffin embedding, and sectioning at 8 µm. Masson's Trichrome staining was performed according to the manufacturer's protocol (Richard Allen).

Limbs were fixed in 10% neutral buffered formalin overnight at 4°C and washed three times in PBS for 10 min. Limbs were then placed on a rocker overnight in 30% acetic acid/70% ethanol/0.3% Alcian Blue stain. When skeletal elements were visibly stained, they were treated with 0.1% trypsin in saturated sodium borate until clear. Some limbs were then treated with Alizarin Red in 1% KOH, then rehydrated in an ethanol series (100%, 95%, 70% and ddH₂O) and run through a 1% KOH/glycerol series of 3:1, 1:1, 1:3 and imaged using a Leica M165 FC stereomicroscope.

Juvenile axolotls 8.8 cm total length (TL) (high=10.1 cm, low=7.4 cm) and 4.58 cm average snout to vent length (SVL) received forelimb amputations at the distal zeugopod. Between days 7 and 14 dpa, individually housed animals were dosed with RA, LE135 or DMSO (n=16/treatment). Drugs were changed every other day. Blastemas containing as little stump tissue as possible were collected from all 48 animals at 14 dpa and single forelimbs from four separate individuals were pooled together to yield four independent biological replicate samples for each treatment group. Total RNA was extracted using the Qiagen RNeasy Kit following the manufacturer's instructions. RNA quality was assessed using an Epoch microplate spectrophotometer, gel electrophoresis, and a 2100 Agilent Bioanalyzer. RNA samples were processed and hybridized to custom A. mexicanum (Amby_002) Affymetrix GeneChips (Huggins et al., 2012) at the University of Kentucky Microarray Core. Expression values were generated using the Robust Microarray Average (RMA) algorithm (Irizarry et al., 2003) and data analysis was performed using the limma software package (Ritchie et al., 2015) in the R environment, generating overall significance statistical values and pairwise comparisons between groups. Venn diagrams were generated using significance values generated for RA/DMSO and LE135/DMSO using the VennDiagram package (Chen and Boutros, 2011). Hierarchical clustering was performed on all 327 significantly changed genes using Cluster (de Hoon et al., 2004) after Log2 transforming the data and mean-centering. Pearson's correlation and average linkage were used to generate a similarity matrix. Trees were visualized using Java TreeView (Saldanha, 2004).

Real-time quantitative PCR collection times were the same as the microarray and biological replicates included four RA-treated samples, four LE135-treated samples and three DMSO-treated controls. cDNA was generated using the Thermo Verso cDNA Synthesis Kit and qPCR with gene-specific primers was performed with ABI PowerSYBR Green PCR Master Mix on a Step-One Plus system following the manufacturer's recommendations. Primers used were: Cyp26a1_F GTGTACCCCGTGGACAATCT, Cyp26a1_R TGCTATGGGTGTTGGGTTTA; Cyp26b1_F CCCTGCTGTAATGGAAGGAT, Cyp26b1_R CGAAGGGCACAATAGGTTTT; Aldh1a1_F AAGACATCGACAAGGCACTG, Aldh1a1_R CCAAAAGGACACTGTGAGGA; Aldh1a2_F GCCAAGACGGTCACAATAAA; Aldh1a2_R CATTCCTGAGTGCTGTTGCT; RARA_F ATACTTGGCAGCCAGAAGGT, RARA_R GCCAACGTTGTATGCATCTC; RARB_F AAAACTCTGAGGGGCTTGAA, RARB_R CTGGTGTGGATTCTCCTGTG; RARG_F CTTCTGCGTTTGATCCTTCA, RARG_R AGTGAGTATGGGGCTGTTCC. Genes were normalized to the control gene FCGBP, which was selected as unchanged in the microarray experiment (primers: FCGBP_F GTTTATGTGGCAGCCTCTCA, FCGBP_R GCCAGCATTAGCTGTGATGT). ΔΔCt was used to calculate fold changes from DMSO controls using the average ΔCt value for each sample.

Treated and control forearms (n=4) were skinned, fixed for 24 h in 10% buffered formalin and then incubated for 24 h in 70% ethanol at room temperature. They were then stained in a 1% phosphotungstic acid/70% ethanol solution for 24 h. The limbs were scanned in the same solution using a microcomputed tomography system (μCT 35, Scanco Medical) (Doube et al., 2010). Scans were acquired with an isotropic resolution of 6 μm, an integration time of 400 ms and a power of 55 kVp. Using BoneJ, we determined the length and the cross-sectional area at midshaft for the radius and the ulna. We reconstructed 3D images of the radius and ulna with the software Mimics.

We thank the many undergraduate volunteers in the Monaghan lab for animal care and discussions about the study.