Conservation analysis of core cell cycle regulators and their transcriptional behavior during limb regeneration in Ambystoma mexicanum

has been one of the major experimental models for the study of regeneration during the past 100 years. Axolotl limb regeneration takes place through a multi-stage and complex developmental process called epimorphosis that involves diverse events of cell reprogramming. Such events start with dedifferentiation of somatic cells and the proliferation of quiescent stem cells to generate a population of proliferative cells called blastema. Once the blastema reaches a mature stage, cells undergo progressive differentiation into the diverse cell lineages that will form the new limb. Such pivotal cell reprogramming phenomena depend on the fine-tuned regulation of the cell cycle in each regeneration stage, where cell populations display specific proliferative capacities and differentiation status. The axolotl genome has been fully sequenced and released recently, and diverse RNA-seq approaches have also been generated, enabling the identification and conservatory analysis of core cell cycle regulators in this species. We report here our results from such analyses and present the transcriptional behavior of key regulatory factors during axolotl limb regeneration. We also found conserved protein interactions between axolotl Cyclin Dependent Kinases 2, 4 and 6 and Cyclins type D and E. Canonical CYC-CDK interactions that play major roles in modulating cell cycle progression in eukaryotes.


Introduction
In multicellular eukaryotes, most developmental processes rely on the transition between the proliferative capacity and the differentiated status of groups of cells.Processes such as embryonic development, postembryonic organ and limb patterning, metamorphosis and regeneration depend on the capacity of cells to divide and to acquire a specific fate.In turn, proliferation and differentiation are influenced by the phase of the cell cycle in which cells are residing.Regeneration of tissues, organs and limbs implies the replacement of populations of lost cells, which requires that specific remnant cells enter into a highly proliferative status.Therefore, the understanding of a given regeneration process must take into consideration the mechanisms by which the cell cycle is modulated in cells that are playing an active role in the multi-stage and complex phenomenon of regeneration (Lange and Calegari, 2010;Filipczyk et al., 2007;Calegari and Huttner, 2003).
Ambystoma mexicanum is an urodele amphibian that has a great capacity for regeneration since it can regenerate lost limbs and organs such as legs, gills, tail, retina, spinal cord, regions of heart and brain (Mitashov, 1997;Cano-Martínez et al., 2010;Simon and Tanaka, 2013).Axolotl limb regeneration is essentially divided into four major stages: wound healing, dedifferentiation, blastema formation and redifferentiation (Yocoyama, 2008).The blastema formation is the main characteristic of epimorphic regeneration and involves the establishment of a heterogeneous population of proliferative progenitor cells.Blastema formation implies that tissue-resident stem cells and differentiated mature cells abandon their G0 status to enter and progress in the cell cycle (Sandoval-Guzmán et al., 2014;Johnson et al., 2018).
Pioneer studies in Ambystoma larva limb regeneration characterized the cell cycle dynamics after amputation by monitoring the incorporation of labeled thymidine (Tassava et al., 1987).A more recent study by Johnson et al. (2018) has shown that 3 days post amputation (dpa) of axolotl limb, epidermal cells that cover the stump reinitiate the cell cycle as revealed by EdU-positive (5-ethynyl-2 ′ -deoxyuridine) nuclei, indicating that cells entered the S-phase.Such experiments also demonstrate that EdU signal dramatically increases at blastema stages (7 and 14 dpa).
EdU-positive cells decrease drastically at redifferentiation stages, suggesting that cells exit the cell cycle prior tissue terminal differentiation.
The E2F-pRb pathway lies at the heart of cell cycle regulation in most multicellular eukaryotes.This pathway not only involves Retinoblastoma family proteins and E2F transcription factors, but also regulatory complexes formed by Cyclins (CYCs) and Cyclin Dependent Kinases (CDKs) proteins.The study of this highly conserved pathway in diverse plants and animals has been extensive, and its crucial role in proliferation, differentiation, cancer and apoptosis has been described, as reviewed in Cross et al. (2011) and Harashima et al. (2013).
In animals, the activity of the Retinoblastoma protein (pRb1) is modulated by post-translational modifications.One of the welldescribed pRb1 modifications is the phosphorylation of diverse residues, exerted by distinct CYC-CDK complexes.Phosphorylation status of pRb1 is central in the dynamics of protein-protein interactions between pRb1 and E2F factors in cell cycle phases (Knudsen and Knudsen, 2008).When pRb1 is hypophosphorylated, binding between the pRb1 pocket domain and the transactivating domain of E2F occurs while in a hyperphosphorylated status pRb1 releases E2F.Free E2F1, in turn, positively regulates the transcription of genes involved in cell cycle progression.In animals, CDK4-CYCD and CDK6-CYCD complexes phosphorylate pRb1 at different sites in mid G1.Then, in late G1, CDK2-CYCE complex phosphorylates additional pRb1 sites (Cross et al., 2011).Phosphorylation of pRb1 by the CDK2-CYCA complex occurs from S phase through G2 and by the CDK1(Cdc2)-CYCB complex from G2 towards M phase (Knudsen and Knudsen, 2008).
Phosphoprotein Phosphatases PP1 and PP2A are responsible for pRb1 dephosphorylation in M phase, which induces cell cycle exit to G0, an event that correlates with the expression of genes that reprograms cells to a differentiated status or maintains some stem cells quiescent (Rubin et al., 2001;Cohen, 2002;Vietri et al., 2006;Kurimchak and Graña, 2015).Also, de novo synthesis of unphosphorylated pRb1 occurs in the resulting daughter cells, which ensures their exit from a proliferative status, unless mitogens and some growth factors, which promote the expression/activity of CYC-CDK complexes, drive cells to re-enter G1.Other important regulators of the cell cycle are the CIP/KIP and CDKN2/INK4 protein families that inhibit the activity of CYC-CDKs leading to a hypo-phosphorylated status of pRb1 thereby preventing cell cycle progression (Knudsen and Knudsen, 2008).Cell cycle entry, with consequent local dedifferentiation and active proliferation of cells, is pivotal for the initiation and progress of limb regeneration.However, studies unraveling the role of cell cycle regulators at the molecular level in this process are very scarce.One of the few studies in this context reported that in cultured myotubes from newt limb, pRb1 phosphorylation by CDK4 and CDK6 is required for cell cycle re-entry (Tanaka et al., 1997).
One of the limitations for the study of the molecular networks that modulate cell cycle and differentiation during axolotl limb regeneration has been the lack of reliable genetic datasets.With the advent of nextgeneration sequencing (NGS), it has recently been possible to obtain the full genome sequence and annotation of this species (Nowoshilow et al., 2018).Besides the full genome sequencing, various transcriptomic approaches have generated very useful information to search and analyze putative functional orthologs for several genes of interest in axolotl (Stewart et al., 2013;Bryant et al., 2017;Caballero-Pérez et al., 2018).
To establish a molecular framework for the major regulators of cell cycle in axolotl, we have identified several putative functional homologs of such genes not only in this species, but also in other Ambystoma species.We have also analyzed their conservation at gene and protein levels and profiled the expression patterns of the respective axolotl mRNAs at 10 different time-points during limb regeneration.Furthermore, we demonstrate in vitro a conserved interaction between specific axolotl CDKs and CYCs.We consider that this molecular framework will be useful and fundamental for researchers aiming to understand the role of cell cycle modulation throughout the limb regeneration process in axolotl (Fig. 1).

Conservation of core cell cycle regulatory genes in axolotl
The genomic data for the axolotl has been quite scarce until recent years and the full genome was sequenced and released just recently (Nowoshilow et al., 2018).Also, a myriad of RNAseq approaches have generated very useful information to search and analyze putative functional orthologs of almost any axolotl gene of interest (Stewart et al., 2013;Bryant et al., 2017;Caballero-Pérez et al., 2018).
As a first step towards the axolotl cell cycle molecular framework, we collected the protein sequence of the most relevant core cell cycle regulators from mouse and human databases since in these two animal species the cell cycle has been thoroughly studied.The obtained sequences were then used to identify putative homologous proteins in diverse databases from vertebrate and invertebrate metazoan species.Among invertebrates, we included Hydra vulgaris and Schmidtea mediterranea, two organisms widely used as model systems for the study of regeneration.Among the vertebrates, in addition to Mus musculus and Homo sapiens, we searched in databases available for 4 species of the Ambystoma genus (A. mexicanum, A. velasci, A. andersoni and A. maculatum), Xenopus laevis, Xenopus tropicalis, Anolis carolinensis, Alligator mississippiensis and Danio rerio.
By using diverse, in silico, approaches (see methods for details), we explored the presence of homologs for several key cell cycle regulators across all the metazoan species investigated (Fig. 2A).Most of the core cell cycle regulators are conserved among vertebrates including the Ambystoma spp.Pocket proteins (pRb1 p105 , pRbl1 p107 and pRbl2 p130 ), Cyclins (CYCA1, CYCA2, CYCB1, CYCB2, CYCB3, CYCD1, CYCD2, CYCE1, CYCE2 and CYCH), E2F transcription factors (E2F1 to E2F8), Cyclin-Dependent Kinases (CDK1 Cdc2 , CDK2, CDK4, CDK6 and CDK7), Protein Phosphatases (PPP1CA and PPP2CA), and two families of CDK inhibitors that include p15, p16, p18, p19, p21, p27 and p57.It is important to note that while p16, a CDK inhibitor of the CDKN2/INK4 family is present in mammals and birds, we were not able to identify any homolog for this protein in reptiles, fishes, amphibians, and invertebrates (Fig. 2A).In contrast to vertebrates, we were not able to identify ortholog sequences for pRb1 p105 , pRbl2 p130 , CYCD1, E2F1, E2F3, E2F4, E2F6, E2F7, CDKN2/INK4 and CIP/KIP in the datasets interrogated for hydra and planaria sequences (see methods).It is also important to note that in planaria, a single E2F transcription factor was found (more similar in sequence to human E2F5), while in hydra there is an expansion of the family with three genes coding for putative homologs of E2F2, E2F5 and E2F8.In the case of E2F3, Fig. 2A refers to the E2F3A isoform, while the E2F3B isoform could not be found in any Ambystoma spp.
Most of the members of CDKN2/INK4 (p15 INK4B , p16 INK4A , p18 INK4C and p19 INK4D ) and CIP/KIP (p21 Cip1/Waf1/Sdi1 , p27 Kip1 , and p57 Kip2 ) protein families are absent in planaria and hydra, with the exception of p27 which is absent in planaria, but present already in hydra.The cell cycle inhibitor p16 is a special case, since this protein is absent in the genomes of the Amystoma spp.and other amphibians analyzed.p16 homolog is also absent in the genomes of reptiles, fishes and invertebrate species analyzed (Gilley and Fried, 2001;Kim et al., 2003).With the exception of p16 and E2F3B isoform, which seem to be absent in axolotl and other Ambystoma spp.analyzed, the rest of the proteins searched have homologs with those of human and mouse.
We found that there is a single pRb-like protein encoded in planaria and hydra genomes, while the vertebrate genomes analyzed, including those of Ambystoma spp., encode for three different pRb proteins, which are homologous to human pRb1 p105 , pRbl1 p107 and pRbl2 p130 (Fig. 2A).Indeed we identified all three proteins for X. tropicalis and X. laevis, while evolutionary analyses of Cao et al. (2010) reported only one pRb protein in X. tropicalis.Although the sequence of pRb-like in hydra and planaria seems to be more similar to that of AmpRbl1 (p107), the percentage of identity is low.This is in agreement with previous studies suggesting that ¨modern¨pRb1, pRbl1 and pRbl2 resulted from two events of genome duplications during metazoan evolution from the common ancestral retinoblastoma protein or apRb (Liban et al., 2017).

Domain conservation analysis in E2F-pRb pathway proteins
Proteins of the E2F-pRb pathway contain crucial polypeptide domains mediating DNA-binding and interaction with other proteins.We therefore performed a protein domain analysis for several members of the E2F-pRb pathway to further explore functional conservation (Fig. 2B, Supplementary Fig. 1, see methods).
For most of the sequences of Cyclins, CDKs and E2Fs families analyzed, we did not observe any obvious differences in the major protein domains (data not shown).We found that in the case of the pocket proteins (pRb1, pRbl1 and pRbl2) the major domains, previously characterized in mouse and humans, are conserved in the respective homologs from Ambystoma spp.In all three proteins, the small pocket domain (SPD) composed by the A-box and B-box, as well as the C-terminal region which, together with the SPD, comprises the large pocket domain (Fig. 2B).In the case of pRbl1, the only pocket protein present in hydra and planaria, the C-terminal domain is absent in the planaria homolog (Fig. 2B).The red line in pRb1 comparison represents a domain found only in Ambystoma ssp.(A. mexicanum, A. velasci, A. andersoni and A. maculatum) which is predicted as a putative cytoplasmic domain.
Whether this is a true functional domain that would affect the nucleocytoplasmic distribution of the protein remains to be elucidated in future studies.In pRbl1, the red line represents an extra cyclin-like domain at the N-terminus.All the domains reported were predicted by using the InterPro Software (Jones et al., 2014).
In the case of E2F transcription factors (E2F TF), our presence and absence analyses indicate that the only E2F TF present in the invertebrate species analyzed is more similar to the E2F5 homologs of the vertebrate species (Fig. 2A).Our domain analysis showed, that although the protein varies in length between species, the major functional domains (DNA-binding and CC-MB domains) in E2F5 homologous proteins are conserved, including the predicted homologs from A. mexicanum and A. velasci (Supplementary Fig. 1, left).Similarly, E2F1 homologs from A. mexicanum and A. velasci conserve both DNA-binding and CC-MB domains (Supplementary Fig. 1A, right).
For Ambystoma spp.CYCD1 and CYCD2 homologs, the canonical Cyclin-like domain is present and its location at the N-terminal region of the proteins is conserved when compared to their homologs in all vertebrate and invertebrate species (Supplementary Fig. 1B).Interestingly, our analysis could not detect a conserved C-terminal region in the single planaria CYCD homolog, sequence identity of which suggests it is more similar to CYCD2 from vertebrates (Supplementary Fig. 1B).
Collectively, these results show that the key components of the E2F-pRb pathway are highly conserved at the protein domain level.However, deeper analyses of the few putative differences found should be carried out and experimentally tested in future studies.Also, a detailed Fig. 1. Cell cycle regulatory proteins and their roles in the checkpoints of the cycle.After limb amputation in A. mexicanum wound healing occurs and the apical epithelial cap is formed, resident connective tissue cells start to dedifferentiate and stem cells leave its quiescent state.Both types of cells are reprogrammed (Left-up corner) and proliferate to populate the blastema.In metazoans, such reprogramming of cells to a proliferative state depends on the activation of CYCD1 and CDK4/6 by mitogens which inhibit the action of CDKN2/INK4 proteins, which are CYCD-CDK repressors.This molecular dynamics leads cells to enter the cell cycle.CYCD1-CDK4/6 complexes phosphorylate pRb proteins, inactivating its capacity to form a repressive complex with E2F transcription factors.The action of pRb-free E2F TFs in transcribing diverse genes is pivotal for the progression of the cycle to S-Phase.The complete release of E2F to exert such action, depends also on the progressive phosphorylation of pRb, not only by CYCD-CDK complexes, but also by the CYCE-CDK2 complex.The counterpart of CYCE-CDK2 complex is CIP/KIP proteins, which inactivate the complex by protein-protein interactions.pRb hyperphosphorylated status is mandatory for the progression of the cycle in both G1-S and G2-M.At the end of mitosis, pRb is dephosphorylated by PPP1CA and PPP2CA complexes, and de novo synthesis of un-phospho pRB in an environment where CYC-CDKs are not active or present, maintains the daughter cells in G0 and some may differentiate.However, it seems that the cellular environment of the blastema niche leads the daughter cells to enter several rounds of divisions until a gradual differentiation of these cells start to form the diverse tissues of the new limb (Created with BioR ender.com).
analysis of conservation of the residues that have been shown as targets of posttranslational modifications in the human homologs, remains to be done.

Phylogenetic analysis of the Ambystoma spp. pRb protein family
Our results show that the homologs for all the members of the retinoblastoma pocket protein family, pRb1, pRbl1 and pRbl2, are present in the Ambystoma species analyzed (Fig. 2A).The domain analyses show that the major functional domains characteristic for pRb1, pRbl1 and pRbl2 are present and their location seems to be conserved too, as well as in keeping a minimal identity to be considered as conserved domain (Fig. 2B).We then performed a phylogenetic reconstruction of the pRb pocket protein family including sequences of ancestral pRb (apRb) from five Pre-metazoa, and the three pRb proteins from 304 vertebrate species, using Pre-metazoan sequences from choanoflagellates and Thecamonas trahens (Apusozoa) to root collapsed and un-collapsed trees (Fig. 3, Supplementary Fig. 2, respectively).Both Maximum Likelihood phylogenetic trees show that pRb1, pRbl1 and pRbl2 protein sequences from Ambystoma spp.grouped with their respective homologs.In the case of Ambystoma pRbl1 homologs they grouped closer to those of frogs.Interestingly, they also group close to their mammalian homologs, Fig. 2. Conservation of core cell cycle regulatory genes along metazoans.A. Presence and absence resume of core cell cycle regulatory genes in diverse metazoan clades.Where purple filled square means presence, white square means absence, red square means that, although the sequence was not found in the current data available for the species, it may be present since the phylogenetically nearest species have the protein.The percentage inside the square reflects the identity of the protein of each species with respect to the homolog from A. mexicanum.*It means that the sequence is not complete.B. Domain conservation analyses of the pocket proteins (pRb1, pRbl1 y pRbl2) among invertebrates, and vertebrate species.The A-box (red box) and the B-box (green box) comprehend the small pocket domain.In turn, the small pocket domain together with the C-terminal region (pink box), composing the large pocket domain.The red line in pRb1 represents a predicted putative domain found only in Ambystoma ssp.and such feature is predicted as a putative cytoplasmic domain.In pRbl1, the red line represents an extra cyclin-like domain in the N-terminal.The putative domains were predicted by InterProScan tool (Jones et al., 2014).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) A. Espinal-Centeno et al. even closer than with those of reptiles (Fig. 3).Although Ambystoma spp.pRb1 proteins group closer to the Amphibian pRb1s, they are phylogenetically closer and appear as a sister clade to the pRb1 homolog of Easter Newt (Notophthalmus viridescens) than to those of frogs.Premetazoan species only have a single apRb protein, which is more similar to the pRb1 homologs of the other species analyzed (Fig. 3).This suggests that although pre-metazoan pRb sequences group closer to pRb1 proteins, the three axolotl pRb proteins (AmpRb1, AmpRbl1 and AmpRbl2) resulted from gene duplications events, and divergence, from an ancestral metazoan sequence, as previously reported in other studies, that did not include the Ambystoma spp.sequences (Cao et al., 2010;Liban et al., 2017).

Phylogenetic analysis of the Ambystoma spp. CDKN2/INK4 family
In mammals, members of the CIP/KIP protein family bind CDKs and inhibit the function of CYC-CDK complexes therefore are inhibitors of cell cycle progression (Besson et al., 2008;Fahmi and Ito, 2019).This family is composed of proteins p21 Cip1/Waf1/Sdi1 , p27 Kip1 , and p57 Kip2 (Sherr and Roberts, 1999;Fahmi and Ito, 2019).The other protein family which has been shown to inhibit the proper function of CYC-CDK complexes is the CDKN2/INK4 family, composed in several metazoan species by p15 CDKN2B , p16 CDKN2A , p18 CDKN2C and p19 CDKN2D .Notably, members of this family bind CYCs instead CDKs, to exert their inhibitory function (Sherr and Roberts, 1999;Besson et al., 2008;Fahmi and Ito, 2019).Due to their importance as inhibitors of the cell cycle, we decided to perform a phylogenetic reconstruction for the members of the CDKN2 protein family.Our collapsed and uncollapsed Maximum Likelihood trees (Fig. 4 and Supplementary Fig. 3, respectively) show that p18 proteins from Ambystoma spp. group together next to their homologs from amphibians.In contrast, Ambystoma spp.p19 homologs grouped together, but more closely to their homologs from fishes and mammals than to those of amphibians.Within the CDKN2 protein family tree, p18 and p19 form a monophyletic clade that could be the result of gene duplication in the last common ancestor of chondrichthyes.The fact that we do not see a p19 sequence from chondrichthyes could be evidence of gene loss in this group of organisms or could be explained because of the absence of a p19 sequence in the interrogated databases.Furthermore, p15 and p16 sequences form a monophyletic clade and within it, the putative p15 homologs found in Ambystoma spp.grouped more closely with p15/p16 homologs from mammals.An interesting observation derived from our presence/absence analysis (Fig. 2A) is the absence of available sequences for p16 protein in the invertebrate species analyzed in this study, as well as absent in Ambystoma spp., birds, and reptiles (Fig. 4).Our findings are in agreement with previous studies that have shown the absence of a functional gene coding for p16 in chicken and newt (Kim et al., 2003;Tanaka et al., 1997).

Quantitation of mRNA levels of core cell cycle regulators throughout axolotl limb regeneration
Transcription of genes that regulate the cell cycle is dynamic and the analysis of their fine-tuned regulation is pivotal to understand the cell reprogramming events that take place during limb regeneration.To explore the transcriptional behavior of the major core cell cycle regulators in A. mexicanum limb regeneration, we isolated RNA from tissues collected at 10 different time points during the process.RNAs were isolated for each time point in three biological replicates and cDNA for each time point was synthesized.RT-qPCR analyses were carried out for each cDNA with a technical triplicate (see Methods for details), to determine the changes in expression at each time point in comparison with time 0.
For a more accurate description of the transcriptional behavior of the 29 genes included in this analysis, we grouped them in families.In the case of CYCs-CDKs, considered as cell cycle activators which are necessary to be expressed in specific checkpoints of the cycle.Our results show that the most upregulated gene corresponds to AmCDK2, which transcript levels increase since 3 dpa, but reach a dramatic change in expression at 10 dpa, a regeneration stage that corresponds to medium blastema (Fig. 5A).Indeed, AmCYCD1, AmCYCE, AmCDK4 and AmCDK6 also reach their top levels of expression at 10 dpa, while transcript levels of AmCYCA2, AmCYCD2, AmCDK1 and AmCDK7 are higher at 7 dpa, at early blastema stage.This shows a high correlation of expression among most AmCYCs and AmCDKs, especially at blastema stages.In fact, the homologs of these genes, in other species, are well known to interact and form specific complexes that promote cell cycle progression, several of them, by phosphorylating pRb proteins.
Acting upstream of CYC-CDK complexes are proteins of the CIP/KIP and CDKN2/INK4 families, which are known to be negative regulators of the cell cycle by binding to either Cyclins or CDKs, inactivating their complexes and in turn, modulating the function of pRb proteins (Fig. 5B).Regarding the expression patterns of genes of the CIP/KIP and CDKN2/INK4 families, we observed that the transcript levels of Amp15 and Amp21 increased since 12 hpa and along all stages, reaching Amp15 a peak at 7 dpa.The transcript levels of Amp18, Amp19, Amp21 and Amp57 genes reached peaks of expression at 10 dpa.On the other hand Amp27 expression did not change along regeneration, showing low transcript abundances (Fig. 5B).
We also found that AmPPP1CA and AmPPP2CA transcript levels increase at 5, 10 and 28 dpa, with peaks of expression at 10 dpa.The evaluated transcripts correspond to the axolotl homologs of the catalytic subunits of PP1 and PP2A complexes.In monkeys both PP1 and PP2A complexes have been implicated in pRb1 dephosphorylation.However, several studies have shown that PP1 is the major pRb1 cell cycle-related phosphatase (Ludlow et al., 1990).
The E2F TFs play a central role in the control of cell cycle, since they are major transcriptional regulators of gene expression in the diverse phases of the cycle.It has been shown in mammals and plants that some E2Fs act as transcriptional activators, while others act as repressors.It has also been shown that pRb1 inhibits the activity of E2Fs through protein-protein interactions, forming transcriptional repressor complexes.In this study, we identified 8 genes coding for E2F TFs in axolotl and named them according to their similarity with the human homologs as AmE2F1, AmE2F2, AmE2F3, AmE2F4, AmE2F5, AmE2F6, AmE2F7 and AmE2F8 (Fig. 2A).In humans there are also 8 genes coding for E2F TFs, and the E2F3 locus codes for two isoforms (E2F3A and E2F3B).Previous studies have shown that E2F1, E2F2 and E2F3A are transcriptional activators, while E2F3B, E2F4 and E2F8 act as transcriptional repressors in humans and mice (Cross et al., 2011;Cuitiño et al., 2019).In the RNAseq datasets and protein databases available for Ambystoma ssp.we found a homolog only for E2F3A, but not for E2F3B.Our RT-qPCR approach shows that, with differences in the levels, all AmE2F genes are transcriptionally upregulated at 3 dpa (Fig. 5C), being AmE2F6 and AmE2F7 the most upregulated.The expression pattern of AmE2F1 shows a maximum of expression at 10 dpa, when blastema is formed and growing and also most AmE2F transcript levels increase when compared to time 0 dpa.While AmE2F3 and AmE2F8 show similar expression patterns, their transcript levels increase since 3 dpa and keep expressed at the remaining time points tested (Fig. 5C).
Finally, quantitation of the relative expression of mRNA levels for members of the pRb family shows that AmpRb1, AmpRbl1 and AmpRbl2, albeit with differences in transcript levels, reach their maximum expression at 10 dpa.Surprisingly, AmpRb1 is the most upregulated at 10 dpa, followed by AmpRbl2 (Fig. 5D).
The data reported above represents the most complete analysis of the expression pattern for axolotl core cell cycle genes in a developmental process and the first such study from the aspect of regeneration.The expression patterns and the observed correlations in expression shown in Fig. 5 will be revisited later in the discussion section.

Protein-protein interaction assays reveal potential AmCYCD1-AmCDK and AmCYCE1-AmCDK complexes
In diverse animal models, consecutive phosphorylation of the pRb1 protein by CYCD1-CDK4/6 and CYCE1/CDK2 complexes in G1 phase brings about the release and activation of E2F transcription factors and promotes cell cycle progression.Because of their critical role we thought to investigate whether the respective axolotl proteins are capable of forming such canonical protein kinase complexes.To test the formation Fig. 4. Maximum Likelihood phylogenetic tree of four proteins from the CDKN2/INK4 family in Metazoa.p15 CDKN2B /p16 CDKN2A , p18 CDKN2C and p19 CDKN2D subfamilies are colored in orange, blue, and magenta, respectively.Only bootstrap values >40 are shown on branches.The NF-kappa-B inhibitor alpha (NKBA in the tree) sequence was used as an external group.In bold we show the Ambystoma ssp.proteins.See Supplementary Fig. 3 for an un-collapsed tree.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Maximum Likelihood phylogenetic tree of four proteins from the CDKN2/INK4 family in Metazoa.p15 CDKN2B /p16 CDKN2A , p18 CDKN2C and p19 CDKN2D subfamilies are colored in orange, blue, and magenta, respectively.Only bootstrap values >40 are shown on branches.The NF-kappa-B inhibitor alpha (NKBA in the tree) sequence was used as an external group.In bold we show the Ambystoma ssp.proteins.See Supplementary Fig. 3 for an un-collapsed tree.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) of complexes between AmCYCD1 and AmCDK2, AmCDK4 and AmCDK6, we transiently co-expressed epitope-tagged AmCYCs and AmCDKs in Arabidopsis protoplasts.Co-immunoprecipitation assays (Co-IPs) showed that AmCYCD1 very weakly binds AmCDK2 but strongly binds to AmCDK4.The interaction was also observed between AmCYCD1 and AmCDK6, although the signal was a lot weaker when compared to that with AmCDK4 (Fig. 6A).
A similar experimental approach was applied to test the binding between AmCYCE1 and AmCDK2, AmCDK4 and AmCDK6.We observed that AmCYCE1 strongly binds to AmCDK2.Interestingly, we also observed that AmCYCE1 blinds to AmCDK4 and AmCDK6, although these interactions appeared to be weaker compared to AmCYCE1-AmCDK2 binding (Fig. 6B).
Collectively, our Co-IPs and the transcriptional analyses suggest that AmCYCD1 may form a functional complex preferentially with AmCDK4, while AmCYCE1 possibly forms a functional complex with AmCDK2 similar to human and mouse.

Discussion
The present study represents the most complete catalog of cell cyclerelated genes in Ambystoma mexicanum reported so far.We identified their presence in axolotl and other Ambystoma spp., and analyzed their conservation along diverse animal clades.This was achieved based on detailed analyses of both genomic and transcriptomic datasets.Moreover, our study describes the transcriptional behavior of most of these genes during the axolotl limb regeneration process.Such transcriptional landscapes include the expression analysis of several genes, which were not explored before in the context of axolotl limb regeneration process.Since a fine tuned regulation of the transitions between cell proliferation and differentiation governs the regeneration process, our findings are pivotal to shed light on the molecular mechanisms that modulate the cell cycle in axolotl.
Previous studies have focused on the evolutionary aspects of cell cycle regulators across the animal kingdom and the resulting data has been extensively discussed (Medina et al., 2016;Liban et al., 2017).Here we discuss from the perspective of the putative homologs of core cell cycle regulators analyzed from Ambystoma spp., and their expression patterns in A. mexicanum limb regeneration.
Our results showed that p16 transcripts are not present neither in the RNA-Seq datasets analyzed nor in the reported protein sequences for any of the Ambystoma spp.(and other amphibians), fishes, reptiles and birds analyzed.This is an interesting observation, since it has been shown for some species that the absence of p16 is not due to the lack of a gene coding for it, but to alterations in genomic arrangements and alternative splicing of the resulting transcripts.Regarding this, in the human genome p15 CDKN2B and p16 CDKN2A genes are encoded in the CDKN2B-ARF-CDKN2A locus, which is composed by the genes that code not only for p15 and p16 mRNAs, but also for ARF (ADP ribosylation factor) mRNA.Indeed HsARF mRNA is specified by the second exon of HsCDKNA via an alternative ORF.A similar genomic arrangement and transcription has been reported for mouse p15, p16 and ARF mRNAs (Kim et al., 2003;Szklarczyk et al., 2007).However, Kim et al. (2003), demonstrated that, in chicken, a similar locus does not encode for p16 CDKN2A .Our results are also in agreement with a previously proposed model in which it was hypothesized that p16 might not be produced in the urodele cells (Tanaka et al., 1997).These observations are interesting since it has been shown in mammals that p16 acts upstream the CYCD-CDK-pRb1 pathway by inhibiting the action of CDK4 and CDK6 and leading to a hypophosphorylated, therefore repressive, status of pRb1.Active pRb1 can then promote partial or terminal differentiation of proliferative cells.Under such scenario, the predicted absence of p16 in urodeles, generated in an epoch when almost no genomic data for urodele amphibians was available, led us to speculate that this might be a major difference between mammals with a limited capacity to regenerate, and urodele amphibians, with high regeneration capacity (Haas and Whited, 2017).Although we cannot exclude that the absence of p16 in axolotl, and other Ambystoma spp., may influence the molecular program that potentiates the regeneration capacity in axolotl in comparison to humans or mouse, we speculate that p16 absence alone is not sufficient to explain the high regenerative ability of the axolotl, since p16 is also absent in birds, reptiles and fishes, species that do not possess the extraordinary regenerative capacity of the axolotl and other urodele amphibians.Also, we can not discard that the regulatory function that p16 protein exerts on its targets could be modulated by other players from the CDKN2/INK4 family.Therefore, future and detailed experiments should be done at the mechanistic level to explore these ideas.Another interesting finding, resulted from our presence-absence analysis, is the absence of the E2F3B isoform in axolotl, and in any of the Ambystoma spp.datasets analyzed.However, we were able to Fig. 6.Protein-Protein interactions between AmCDKs and AmCYCD1-E1. A. Interactions between AmCYCD1-AmCDK2, AmCYCD1-AmCDK4 and AmCYCD1-AmCDK6.B. Interactions between AmCYCE1-AmCDK2, AmCYCE1-AmCDK4 and AmCYCE1-AmCDK6.Proteins were expressed in Arabidopsis protoplasts.Complexes immunoprecipitated with anti-c-Myc antibodies were analyzed on protein gel blots using anti-HA and anti-c-Myc antibodies to detect protein-protein binding (upper panels in A and B) and expression of AmCDKs (lower panels in A and B) respectively.identify in all Ambystoma spp.analyzed, the homologous sequences coding for the E2F3A isoform.In mammals, E2F3A acts as a transcriptional activator, while E2F3B has been shown to act as a repressor (Asp et al., 2009).This is a quite interesting observation since it has been shown that E2F3A and E2F3B exert opposing regulation of SOX2 in mouse neural precursors (Julian et al., 2013), and axolotl SOX2 is critical for promoting the proliferation of neural stem cells upon tail amputation (Fei et al., 2014).Weather AmE2F3A acts regulating AmSOX2 transcription and the absence of AmE2F3B is of relevance on AmSOX2-mediated tail regeneration, should be experimentally addressed in the future.
Developmental processes, such as limb regeneration in multicellular organisms, rely on finely regulated cell proliferation which, in turn, critically depends on cell-cycle-dependent mRNA oscillations.It has been proposed that such oscillatory behavior of cell-cycle-regulated genes are partially related to E2F-dependent transcriptional activation and repression (Giangrande et al., 2004;Li et al., 2008;Cuitiño et al., 2019).
Several interesting observations resulted from our expression pattern analyses.The transcriptional analyses showed that AmCDK2, AmCDK4 and AmCDK6 reach their maximum peak of expression in blastema stages, specifically at 10 dpa.Surprisingly, AmCDK2 is the most upregulated gene at 10 dpa, when compared to time 0 dpa.This suggests that several cells in this time point are preparing to enter S-phase (Fig. 1).AmCDK2 transcription is upregulated since 24 hpa, as well as its putative complex partner AmCYCE.Also, AmCDK4 and AmCYCD2 transcript levels increase since early times after amputation (Fig. 5A).The early transcriptional activation of axolotl CYCs and CDKs, since 12 hpa, led us speculate on two probable scenarios: i) That existing proteins after amputation may form axolotl CYC-CDK complexes which phosphorylate AmpRb and AmE2F may induce the further transcription of AmCYCs and AmCDKs in further timepoints or ii) The transcriptional activation of these genes is independent of AmpRb status, for example the transcriptional activation of CYCs could be caused by mitogens such as growth factors that act through the MAP-Kinase signaling pathway, as previously shown in other models and experimental systems (Roovers et al., 1999).Also, our expression analyses showed that all AmE2F tested are transcriptionally induced at 3 dpa (Fig. 5C).By 3 dpa the wound healing process is ending and de-differentiation may start to take place, and cells start to leave G0 status and enter G1.These expression patterns coincide with measurements made in diverse phases of the cell cycle using mouse embryonic fibroblast (MEF) cell culture, where E2F1, E2F2, E2F3A, E2F7 and E2F8 mRNA levels were higher in G1 phase cells when compared to cells in G0, and their expression peaks at S and G2 phases (Cuitiño et al., 2019).AmE2F1 is the most upregulated gene of this TF family, with transcript levels reaching a maximum of expression at 10 dpa.Moreover, after 3 dpa the expression levels of AmE2F3, AmE2F5, AmE2F6 and AmE2F7 remain relatively constant along the process, with exception of another peak of expression of AmE2F6 and AmE2F7 at 21 dpa.Interestingly, mRNA levels of E2F3b, E2F4, E2F5 and E2F6 did not change across all cell-cycle phases in MEF cell cultures (Cuitiño et al., 2019).
AmE2F1 and AmE2F8 homologs, which in human and mouse act as transcriptional activators and repressors, respectively, are the most upregulated E2Fs at blastema stages (7 and 10 dpa).These transcriptional patterns suggest that in a population of cells that constantly proliferate, as in blastema stages, both transcriptional activation and repression occurs.Indeed it has been hypothesized that transcriptional repression mediated by E2F8 in mid-S-G2 allows a temporary and repressive state that is necessary for a next wave of cell-cycle-dependent gene expression to be reactivated during the next round of cell cycle (Cuitiño et al., 2019).It has been shown that E2F8 lacks the ability to form a repressive complex with pRb1, since mouse and human E2F7 and E2F8 lack the DP-dimerization and the retinoblastoma binding domains (Maiti et al., 2005).Whether such mechanisms of cell cycle control are occurring during axolotl limb regeneration needs to be demonstrated via future functional experiments.
Other interesting observations, which resulted from our transcriptional analysis, is the high expression of diverse CYC-CDK repressors.It is noteworthy the transcriptional behavior of Amp21 which reaches its highest expression levels 10 dpa.Also, other CDK repressors such as Amp57 and Amp18 are upregulated mainly at 10 dpa.In contrast, Amp15 expression is upregulated from 24 hpa until 7 dpa, and its expression is down regulated at 10 dpa.This indicates that cell cycle progression is dynamically modulated, suggesting that CIP/KIP and CDKN2/INK4 proteins are being produced to inhibit CDK2, CKD4, CDK6 and preparing cells to exit the cycle and perhaps, help in the differentiation of certain cells in specific stages of the regeneration process.
It should also be taken into consideration that transcript levels may not coincide with levels of the respective encoded proteins and their activity.Several cell cycle regulators have been shown to interact with each other, and such interactions are dependent on posttranslational modifications.This is true for pRb proteins, CYCs, CDKs and E2Fs.Whether critical modifications are conserved in the putative functional homologs in axolotl and whether modifications are part of the mechanisms that govern cell reprogramming events during regeneration, need to be experimentally explored in future studies.
As a proof of concept for this, we explored the interaction between axolotl CYCs and CDKs.CYCD1-CDK2/4/6 and CYCE1-CDK2/4/6 complexes are two of the most conserved protein-protein interactions during plant and animal evolution.The specificity of the formed complexes and their pivotal roles in diverse stages of cell cycle have been extensively described as major switches for cell proliferation, therefore underlying a myriad of developmental processed and morphogenetic events (Cao et al., 2010;Zhu and Pearson, 2013;Medina et al., 2016).Our Co-IP assay using Arabidopsis protoplast system has shown very clean results on protein-protein interaction, this system is quite clean and avoids the use of axolotl cells, however we should consider that the proper interactions in vivo in the axolotl remain to be confirmed.Future experiments are needed in order to show if the strength and specificity of these interactions is similar in axolotl tissues, where several other axolotl cell cycle players might be influencing these CYCD1-CDK interactions.
Our results showed that AmCYCD1 strongly binds to AmCDK4 (Fig. 6A).Interestingly, the transcriptional patterns of AmCYCD1 and AmCDK4 suggest co-expression, since they grouped in the Heatmap visualization of our RT-qPCR analyses (Fig. 5A).We also observed that AmCYCE1 preferentially binds to AmCDK2 (Fig. 6B).A weak interaction was also observed between AmCYCE1 and AmCDK4, and AmCDK6.The transcriptional patterns of AmCYCE and AmCDK2 grouped in the Heatmap visualization of our RT-qPCR analyses (Fig. 5A).Moreover, the relative expression levels of AmCDK6 transcripts are low throughout regeneration, and less abundant than those of other AmCDK genes, suggesting that the most probable AmCYCE partner in vivo is AmCDK2.
These findings reveal highly conserved interactions, since in human and mouse CYCD1 forms complexes with CDK4 and CDK6, to phosphorylate pRb at the G1-S checkpoint, while CYCE partner is CDK2 and work together to continue pRb phosphorylation in late G1-S, which is a major event that brakes the interaction between pRb and E2Fs TFs, allowing E2Fs to exert their role on gene transcription activation for proper cell cycle progression (Cao et al., 2010;Cuitiño et al., 2019).
pRb proteins are master regulators of the cell cycle in diverse kingdoms of life and their role strongly depends on their phosphorylation status which, in turn, relies on the action of CYC-CDK complexes.In this study, genes coding for putative functional homologs of the three pRb proteins (pRb1, pRbl1 and pRbl2) were identified in 4 Ambystoma spp., including the axolotl.Phylogenetic and protein domain analyses revealed conservation with their homologs in the amphibian clade.Our phylogenetic analyses show that Ambystoma spp.pRb1 sequences group closer to the pre-metazoan sequences (apRb).Ambystoma pRb1, pRbl1 and pRbl2 might follow a similar evolutionary path than that described in previous studies, which indicate that these three proteins resulted A. Espinal-Centeno et al. from gene duplications and divergence of an ancestral metazoan sequence (Cao et al., 2010;Liban et al., 2017).
One of the model organisms with the highest regeneration capacity in the animal kingdom is Schmidtea mediterranea.The single pRb protein in fresh-water planaria, Smed-Rb, is more similar to p107 homologs and is highly expressed in proliferative cells.The characterization of Smed-Rb RNAi lines has shown that it is required for cell-cycle progression, while it was not essential for proper differentiation (Zhu and Pearson, 2013).This correlation between the absence of a canonical pRb1 protein in planaria and its high regeneration capacity led us to speculate that there might be a partial loss of the role of AmpRb1 as a canonical repressor of the cycle, which together to the absence of p16 and E2F3B homologs, may be related with the high regenerative capacity in axolotl.To test such hypotheses more detailed conservation analyses on axolotl pRb proteins, at the level of motifs and residues, should be done to reveal the existence of subtle differences among pRb proteins from axolotl and their homologs in species with a lower regeneration capacity.
The changes in transcript levels, observed in Fig. 5, are subtle for several of the genes tested, and in some cases seem contradictory with the processes occurring in a regeneration stage.This suggests that some of these genes might not be regulated at the transcriptional level, but protein modifications may be occurring.We should also consider that the change in the expression of a gene might be restricted to a specific subset of cells or cell lineage, and such changes might be masked by the non-changing transcript levels of the same gene in the rest of the tissue used for RNA isolation.It would be very useful to determine the in situ expression of transcripts, and proteins, encoded in these genes in order to generate transgenic markers for cell cycle stages.These novel tools could make possible single cell expression analyses of such cell cycle regulators.
Although our study describes the most complete catalog of cell cycle related genes in axolotl reported so far, as well as a detailed analysis of their transcriptional behavior during limb regeneration, the functional validation for each of their, putatively homologs, genes remains pendant.Therefore, future research must include not only loss and gain of function approaches to evaluate each of the genetic networks and pathways that control cell cycle status, not only during limb regeneration but also during development.For example, the relevance of the absence of p16 and E2F3B, in the context of regeneration, remains to be explored.
Finally, we expect that making public the data described in this study could be of use and support for future functional studies aimed to reveal and understand the molecular mechanisms, related to cell cycle control, underlying axolotl regenerative capacity.The differences in such molecular mechanisms between axolotl and humans are of key relevance for regenerative medicine.

Domain conservation analyses
The orthologous sequence for all 35 cell cycle proteins was scanned in the InterPro program to find domain old and new putative.InterPro uses predictive models for various different databases (Jones et al., 2014).

Animal limb amputation and tissue processing
Juvenile A. mexicanum specimens from Xochimilco, were obtained A. Espinal-Centeno et al. from the Unit of Environmental Management at the Centro de Investigaciones Acuáticas de Cuemanco, Universidad Autónoma Metropolitana of México City, Campus Xochimilco.Three individuals 7.5-8.5 cm size, in length head-tail for each experiment time of regeneration were used.The animals were anesthetized by immersion in a 0.01% solution of benzocaine (ethyl 4-aminobenzoate, SIGMA), using a scalpel the anterior-right limb was amputated at the mid-forearm level.The limb, or regenerated tissues, of each individual were processed for RNA extraction and the obtained RNA was used for cDNA synthesis.
The animals experiments were conducted to fully agree with the Mexican Official Norm (NOM-062-ZOO-1999) "Technical Specifications for the Care and Use of Laboratory Animals" based on the Guide for the Care and Use of Laboratory Animals "The Guide", 2011, NRC, USA, with the Federal Register Number # BOO.02.08.01.01.0095/2014, awarded by the National Health Service, Food Safety and Quality (SENASICA-SADER).The Institutional Animal Care and Use Committee (IACUC) from the CINVESTAV approved the project "Management and husbandry of Ambystoma spp and experimental processing of tissue for functional analyses and genetic expression" ID animal use protocol number: 0209-16.

RT-qPCR quantification
Total RNA samples of different times (0 hpa, 12 hpa, 24 hpa, 3 dpa, 5 dpa, 7dpa, 10 dpa, 14 dpa, 21 dpa and 28 dpa) from three individuals of 7.5-8.5 cm in length, head-tail, for each regeneration time point, were used.Each organism was anesthetized by immersion in a 0.01% benzocaine (ethyl 4-aminobenzoate, SIGMA) solution.Then the anterior right limb was amputated at the mid-forearm level.The tissue of each of the tree organisms were collected in each time point and processed for RNA extraction using TRIzol Reagent (Invitrogen).RNA was used for the synthesis of cDNA according to the manufacturer protocol SuperScript III RNA to cDNA kit (Invitrogen), we used 3 μg of total RNA per 20 μL reaction.The expression levels of specific mRNAs were determined by RT-qPCR using SYBR Green PCR master mix in 20 μL reaction.We used three biological replicates with three independent RT-qPCR reactions for each data point and we included the primers efficiency in the calculations.The primer sequences used in these experiments are in Supplementary Table II.The resulting data for each time point and transcript was normalized with ODC1 expression levels previously reported to axolotl limb regeneration (Guelke et al., 2015).The relative expression was calculated by 2^(− Δ Δ Cт) method, subtracting the expression levels of time 0 to those each time point analyzed.

Fig. 3 .
Fig.3.Maximum Likelihood phylogenetic reconstruction of the retinoblastoma pocket protein family in Metazoa.pRb1 p105 , pRbl1 p107 , and pRbl2 p130 groups are colored in magenta, orange, and blue, respectively.Only bootstrap values >45 are shown on branches.Pre-metazoan sequences from choanoflagellates and Thecamonas trahens (Apusozoa) were used to root the tree.See Supplementary Fig.2for an un-collapsed tree.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Maximum Likelihood phylogenetic reconstruction of the retinoblastoma pocket protein family in Metazoa.pRb1 p105 , pRbl1 p107 , and pRbl2 p130 groups are colored in magenta, orange, and blue, respectively.Only bootstrap values >45 are shown on branches.Pre-metazoan sequences from choanoflagellates and Thecamonas trahens (Apusozoa) were used to root the tree.See Supplementary Fig.2for an un-collapsed tree.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Heatmap visualization of RT-qPCR analysis of mRNA levels of core cell cycle regulators during axolotl limb regeneration.RT-qPCR was performed on cDNA using gene-specific primers for diverse axolotl homologous genes to those involved in cell cycle regulation, along 10 time-points after limb amputation.Relative quantification of each gene expression level was normalized according to the AmODC1 gene expression.A. AmCYCs and AmCDKs B. AmCIP/KIP, AmCDKN2/INK4, and AmPhosphatase, families.C. AmE2F Transcription factors.D. AmpRb1 and AmpRb-like genes.The values were validated for three biological replicates by RT-qPCR, and expressed as relative quantification [2^(-ΔΔCт)], data are the mean of three biological replicates, each one with three technical replicates.Purple represents up-regulation, white represents no change, and yellow represents down-regulation, with respect to time 0. Graduation in colour intensity represents the variations in mRNA expression values, Standard deviations (+/− ) were calculated and presented below the expression value.Limb regeneration stages are described as follows, WH: Wound healing, DD: Dedifferentiation, EB: Early blastema, MB: Medium blastema, LB: Late blastema, ED: Early digit (formation), MD: Medium digit (formation).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)