The Dynamic Landscapes of Circular RNAs in Axolotl,a Regenerative Medicine Model,with Implications for Early Phase of Limb Regeneration

Abstract

Circular RNAs (circRNAs) are of relevance to regenerative medicine and play crucial roles in post-transcriptional and translational regulation of biological processes. circRNAs are a class of RNA molecules that are formed through a unique splicing process, resulting in a covalently closed-loop structure. Recent advancements in RNA sequencing technologies and specialized computational tools have facilitated the identification and functional characterization of circRNAs. These molecules are known to exhibit stability, developmental regulation, and specific expression patterns in different tissues and cell types across various organisms. However, our understanding of circRNA expression and putative function in model organisms for regeneration is limited. In this context, this study reports, for the first time, on the repertoire of circRNAs in axolotl, a widely used model organism for regeneration. We generated RNA-seq data from intact limb, wound, and blastema tissues of axolotl during limb regeneration. The analysis revealed the presence of 35,956 putative axolotl circRNAs, among which 5331 unique circRNAs exhibited orthology with human circRNAs. In silico data analysis underlined the potential roles of axolotl circRNAs in cell cycle, cell death, and cell senescence-related pathways during limb regeneration, suggesting the participation of circRNAs in regulation of diverse functions pertinent to regenerative medicine. These new observations help advance our understanding of the dynamic landscape of axolotl circRNAs, and by extension, inform future regenerative medicine research and innovation that harness this model organism.

Introduction

Circular RNAs (CircRNAs) are stable, conserved across species, abundantly expressed, tissue- and developmental stage specific, and one of the recently discovered RNA types. They were initially considered nonfunctional aberrant RNAs that arise from errors during the gene transcription process. However, studies have shown that circRNAs are essential to control many biological processes by acting as miRNA sponges, gene expression regulators, and interactors of proteins. Regarding their origin, circRNAs are divided into three main groups: exonic circRNAs, circular intronic RNAs, and exonic–intronic circular RNAs (Memczak et al., 2013). They lack the 5′ cap or 3′ poly-A tail and are circularized by covalent bonds (Li et al., 2018).

circRNAs play diverse biological roles in homeostasis and disease pathology through several mechanisms, mainly competing with host genes for splicing machinery, serving miRNA-binding sites, and modulating gene expression at the post-transcriptional level (Kulcheski et al., 2016). For example, ciRS-7 circRNA, also called CDR1as, functions as a typical miRNA sponge and has more than 70 binding sites for miR-7 (Hansen et al., 2013a). The inability of miR-7 to regulate its target mRNAs due to binding to ciRS-7 is observed in many cancers such as lung cancer and astrocytoma (Hansen et al., 2013b; Weng et al., 2017). In addition to being miRNA sponges, circRNAs can act as sponges for RNA-binding proteins (RBPs) to adjust protein levels. Considering the roles that RBPs take in almost all cellular processes, including proliferation, apoptosis, and differentiation, circRNAs have the potential to control almost all biological processes. For example, circMBL, a conserved among invertebrates and vertebrates, can bind to the Mannan binding lectin (MBL) protein from many sites to regulate the various pathways in muscle and nervous tissues (Ashwal-Fluss et al., 2014; Pamudurti et al., 2022).

As another example, the roles of circFOXO3 in the cell cycle and apoptosis are well described. It has been documented that the circRNA of circFOXO3 interacts with p21 and cyclin-dependent kinase 2 to form a ternary complex, which decreases the activity of cell cycle essential genes, Cyclin A and Cyclin E, and thus blocks the cell cycle progression (Du et al., 2016).

The association of aberrant circRNA expression with many diseases further supports their crucial roles in maintaining homeostasis and health. However, circRNA's functions in regeneration, particularly limb regeneration, are under investigated. Salamander axolotl (Ambystoma mexicanum) is an astonishing model organism to explore the molecular basis of organ and extremity renewal due to its exceptional functional restoration capacity (Bölük et al., 2022). Cutting-edge next-generation sequencing technologies have been adopted to dissect the alteration in gene expression profile throughout the regeneration process. Axolotl reference datasets generated from RNA sequencing, proteomics, microbiome, and metabolome analyses have expanded the current understanding of the regeneration concept (Demircan et al., 2020; Demircan et al., 2019; Demircan et al., 2017; Gerber et al., 2018; King and Yin, 2016; Liu et al., 2022; Qin et al., 2021; Rao et al., 2014; Sibai et al., 2019; Wei et al., 2022). Yet, the circRNA data for axolotl are still unavailable. To identify the axolotl's circRNAs and predict their roles in limb regeneration, the circRNA profile during the early phase of limb regeneration has been mapped.

Differentially expressed (DE) circRNAs at 1- and 14-day postamputation (DPA) compared with intact limb were listed to enrich gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome pathways. The findings will contribute toward an extension of circRNA roles and better utilization of axolotl in regeneration studies.

Materials and Methods

Experimental design and animal experiments

The Istanbul Medipol University local Ethics Committee approved all performed animal experiments with authorization number 38828770-E15956. A total of 27 randomly selected sibling wild-type axolotls (1-year-old, 12 cm in length) were used in this study to form 0-, 1-, and 14-DPA time points. Axolotls were obtained from the Istanbul Medipol University Medical Research Center and bred as described elsewhere (Demircan et al., 2017). To have repeatable accuracy, nine animals of each time point group were subgrouped into three biological replicates (R1, R2, and R3). Before amputations and sample collection, animals were anesthetized using 0.1% ethyl 3-aminobenzoate methanesulfonate (MS-222; Sigma-Aldrich, St Louis, MO, USA). Following the amputation, samples were collected at 0-, 1-, and 14-DPA around the cut site and frozen immediately to store at −80°C until RNA isolation.

RNA isolation and circRNA sequencing

TRIzol reagent (Cat. No. 15596018; Invitrogen™) was used to isolate the total RNA from frozen tissues collected at 0-, 1-, and 14-DPA according to the manufacturer's procedure. Previously established protocols (Panda et al., 2017; Xiao and Wilusz, 2019) were followed to purify circRNA from total RNA by depleting the mRNA and rRNAs. Briefly, 20 μg of total RNA samples were mixed with 20 U RiboLock RNase inhibitor (Thermo Fisher Scientific), 1 × RNase R buffer, and 20 U RNase R (RNR07250; Epicenter) to form a 20 μL reaction volume, which was incubated at 37°C for 30 min to eliminate the linear RNA. Following the phenol:chloroform sodium acetate precipitation, a poly(A) tail was added to the residual linear RNAs using a commercially available Polyadenylation Kit (AM1350; Thermo Fisher Scientific) according to the producer's protocol. The Poly(A) Purist™ MAG Kit (AM1922) was used to select the RNAs with a poly(A) tail to discard any remaining linear RNAs.

Since there is no free hydroxyl group at the 3′ end of circRNAs, polyA tails cannot be added to them, and they were not affected by this step. Purifying circRNAs was followed by preparing a sequencing library using the TruSeq Small RNA Library Prep Kit protocol according to the manufacturer's instructions (Cat No: RS-200-0012; Illumina). miRNA libraries were quantified by quantitative polymerase chain reaction (qPCR) using the KAPA Biosystems Library Quantification Kit for Illumina platforms (Cat No: KK4844; Roche). Sequencing libraries were loaded at 20 pM concentration on the Illumina NovaSeq sequencer.

Raw data processing, mapping, and quantification of circRNAs

The Trinity algorithm (Grabherr et al., 2011) was employed to perform de novo assembly of raw reads, resulting in a database of putative axolotl circRNAs. The annotations for these putative circRNAs are computed based on the best hits (lowest e-value) from BLAST+ (Camacho et al., 2009) done against 140,538 human circRNAs from Circbank (Liu et al., 2019) and CircAtlas (Wu et al., 2020) databases. The RNA-Seq by Expectation Maximization (RSEM) protocol (Li and Dewey, 2011) is implemented to compute cirRNA expression profiles based on this putative axolotl circRNA database. First, circRNA expression count vectors are computed for biological replicates at 0-, 1-, and 14-DPA time points. Then, EBSeq (Leng et al., 2013) is used to identify circRNAs DE between 0-DPA versus 1-DPA and 0-DPA versus 14-DPA time points.

Validation of circRNA sequencing results

To validate the obtained DE circRNA list by following a previously published protocol (Vromman et al., 2021), five DE circRNA were randomly selected among 0-, 1-, and 14-DPA. cDNA of rRNA and mRNA-depleted RNA samples used in library preparation were obtained employing random hexamers and the ProtoScript First-Strand cDNA Synthesis Kit (E6300S; NEB) according to the manufacturer's procedure. qPCR was carried out utilizing circRNA-specific primers (Supplementary Table S1) and the SensiFAST SYBR^® No-ROX Kit (BIO98005; Bioline), following the manufacturer's suggestions. The quantitative reverse transcription-PCR (qRT-PCR) of three biological and three technical replicates was conducted using Roche LightCycler, and relative circRNA expression levels were calculated with the 2^−ΔΔCt method (Livak and Schmittgen, 2001).

circRNA target prediction, enrichment analyses, and visualization

“Enhanced volcano” R package was used to display the DE circRNAs with p < 0.01 and |FC| > 2.0 on volcano plots. Lists of annotated circRNA and miRNA–circRNA interactions were downloaded from the Circbank database (Liu et al., 2019), the human circRNA database, and these two lists were merged to get a final data frame of circRNAs and corresponding miRNAs. Then, axolotl circRNA lists were matched with the merged list obtained in the previous step using the common column “CircbaseID” to identify the miRNAs potentially sponged by axolotl circRNAs. The yielded identified miRNA lists were used to find validated mRNA targets of those miRNAs by searching in “mirecords,” “mirtarbase,” and “tarbase” databases using the “multiMiR” R package (Ru et al., 2014). Found genes with EntrezID were used for GO and KEGG enrichment analyses utilizing “clusterProfiler” package in R (Yu et al., 2012). To enrich the Reactome pathways by DE circRNAs, the ReactomePA package (Yu and He, 2016) was employed. In all enrichments, to decrease the false discovery rate, p-values were adjusted with Benjamini and Hochberg procedure, and cutoffs of adjusted p-value and q-value were set to 0.05.

The “org.Hs.eg.db” (human) database was selected as the background gene list. The top enriched pathways and processes were plotted using “ggplot2” and “ggpubr” R packages.

Protein–protein interaction analyses

STRINGapp (based on STRING version 11.5) in Cytoscape was used to construct and visualize the protein–protein interaction (PPI) networks (Szklarczyk et al., 2019). To analyze the clustering and interaction of DE proteins, upregulated and downregulated proteins in 1-DPA and 14-DPA samples were queried. The top 20 DE proteins with the highest degree were picked, and the PPI network of these 20 proteins was drawn.

Results

The raw reads (221,112,171 read pairs in total, M = 24,568,019 SD = 6,476,128) generated in this study were submitted to the NCBI BioProject and SRA databases (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA852678/). The de novo assembly of raw reads resulted in the identification of 35,956 putative axolotl circRNAs, of which 34,811 had at least one orthologous human circRNA based on BLAST+ searches done against Circbank and CircAtlas databases. Principal component analysis (PCA) of the circRNA expression data (Fig. 1A) illustrated time point separation and clustering of replicates within the same group (PC1 = 55.6%, PC2 = 21.1%). PC1 demonstrated a clear separation between 14-DPA and 0-, and 1-DPA, while PC2 seemed to separate 0-DPA and 1-DPA groups. Five thousand three hundred thirty-three unique circRNAs with a human orthology identified in this study distributed across the time points as 3686 circRNA for 0-DPA, 3533 for 1-DPA, and 3286 for 14-DPA. Thirty-seven percent of these circRNAs were commonly found in all time points, whereas 23% were identified in two of the three time points (Fig. 1B).

FIG. 1.

Cluster analysis was performed using circRNA expression data. PCA was generated for three time points (0-DPA, 1-DPA, and 14-DPA) and a timewise separation was observed (A). Venn diagram was generated to provide a visual representation of the similarities and uniqueness of identified circRNAs among samples (B). Volcano plot depicting significantly DE circRNAs (p < 0.01, |FC| > 2) in 1-DPA (C) and 14-DPA (D) samples. Venn diagram was generated to illustrate the number of DE circRNAs in timewise comparisons (E). Expression levels of randomly selected circRNAs. Among the DE circRNAs, five of them were selected to validate circRNA sequencing (F). circRNA, circular RNA; DE, differentially expressed; DPA, day postamputation; FC, fold change; PCA, principal component analysis.

To analyze the RSEM data displaying log fold change against p-value, volcano plots were generated for 0- and 1-DPA, and 0- and 14-DPA comparisons (Fig. 1D, E). In the 0-DPA versus 1-DPA group, a total of 20,197 circRNAs were assessed for significance, and 323 circRNAs met the cutoff criteria (Fig. 1C, E). In contrast, for the 0-DPA versus 14-DPA comparison, out of 19,284 circRNAs analyzed, 606 circRNAs showed significance (Fig. 1D, E). One hundred forty-three circRNAs were shared between two comparisons (Fig. 1E). Validation of the sequencing and in silico analysis was performed using qRT-PCR, and the expression levels of five randomly selected DE circRNAs (hsa-circ-0007197, hsa-circ-0022980, hsa-circ-0068343, hsa-circ-0124257, and hsa-circ-0131531) in the experimental groups were measured and visualized (Fig. 1F). A similar trend between RNA seq and qRT-PCR results stimulated us to continue with downstream analysis.

Next, DE circRNAs without hits in the human circRNA database were filtered, and the miRNAs potentially interacting with the annotated DE circRNAs were investigated (Supplementary Table S2). Among the identified miRNAs, those commonly found in either upregulated or downregulated lists were excluded, and the remaining miRNAs were presented in Supplementary Table S3. At 1-DPA, 23 miRNAs were upregulated, and 79 miRNAs were downregulated, while at 14-DPA, 11 miRNAs were upregulated and 15 miRNAs showed potential downregulation as interactors of annotated DE circRNAs. Validated targets of these miRNAs (Supplementary Table S4 and Fig. 2A) were then analyzed to enrich the biological process (BP) terms, Reactome pathways, and KEGG pathways. In the comparison between 1-DPA and 0-DPA, the miRNAs upregulated and downregulated targeted 231 and 2637 genes, respectively, while in the comparison between 14-DPA and 0-DPA, the upregulated miRNAs targeted 7967 genes, and the downregulated miRNAs targeted 565 genes.

FIG. 2.

The number of genes exhibiting significant upregulation or downregulation in their expression levels on days 1-DPA and 14-DPA were depicted using a Venn diagram (A). For the visualization of the outcomes from the GO enrichment analysis, a dot plot was generated to present the top 10 activated or suppressed BP terms at 1-DPA or 14-DPA (B). BP, biological process; GO, gene ontology.

The top 10 enriched BP terms in 1-DPA samples included “cell cycle,” “nervous system development,” and “nitrogen compound transport,” while the top 3 suppressed pathways were “tRNA metabolic process,” “mitochondrial respiratory chain complex assembly,” and “tRNA processing” (Fig. 2B). In contrast, the top activated BP terms in 14-DPA samples were “cellular response to stress,” “cell cycle,” and “mitotic cell cycle,” while the top suppressed terms were “negative regulation of nitrogen compound metabolic process,” “protein-containing complex subunit organization,” and “protein-containing complex assembly” (Fig. 2B). The altered BP terms are listed in Supplementary Table S5.

Furthermore, the genes enriching KEGG pathways were analyzed (Fig. 3A). The important activated pathways in the 1-DPA samples included “TGF-beta signaling pathway,” “Mineral absorption,” and “Fc gamma R-mediated phagocytosis.” On the other hand, the top 3 suppressed pathways were “Proteoglycans in cancer,” “FoxO signaling pathway,” and “Human T cell leukemia virus 1 infection.” In the 14-DPA group, the enriched pathways included “Proteoglycans in cancer,” “mammalian target of rapamycin (mTOR) signaling pathway,” “MAPK signaling pathway,” and “Autophagy—animal,” while the suppressed pathways were “Apoptosis,” “NF-kappa B signaling pathway,” and “Herpes simplex virus 1 infection.” The lists of upregulated and downregulated KEGG pathways are presented in Supplementary Table S6.

FIG. 3.

To illustrate the top 10 KEGG pathways that exhibited activation or suppression at 1-DPA or 14-DPA following an additional enrichment analysis a barplot was employed (A). The top 10 Reactome Pathways that were significantly activated or suppressed at 1-DPA or 14-DPA were visually represented by Utilizing a barplot (B). KEGG, Kyoto Encyclopedia of Genes and Genomes.

The top 10 activated and suppressed Reactome pathways (Supplementary Table S7) at 1-DPA and 14-DPA are illustrated in Figure 3B. Notable activated Reactome pathways in 1-DPA included telomere length-related pathways, such as “Telomere C-strand (Lagging Strand) Synthesis,” “Extension of Telomeres,” and “Processive synthesis on the C-strand of the telomere,” as well as DNA repair or cell cycle-related pathways like “HDR through Homologous Recombination (HRR),” “Recruitment of NuMA to mitotic centrosomes,” and “Chromosome Maintenance.” Notch signaling-related pathways, such as “NOTCH4 Activation and Transmission of Signal to the Nucleus,” “Signaling by NOTCH4,” and “NOTCH2 Activation and Transmission of Signal to the Nucleus,” were suppressed at 1-DPA. At 14-DPA, the top activated Reactome pathways included “Mitotic Anaphase,” “Mitotic Metaphase and Anaphase,” “Semaphorin interactions,” and “Separation of Sister Chromatids.” Noteworthy Reactome pathways described at 14-DPA were “Cellular Senescence,” “DNA Damage/Telomere Stress-Induced Senescence,” “Diseases of programmed cell death,” and “Senescence-Associated Secretory Phenotype (SASP).”

Finally, the PPI network was examined for the DE proteins identified in the 1-DPA and 14-DPA samples. The top 20 DE proteins with the highest number of interactions are depicted in Figure 4. In the 1-DPA network, upregulated DE proteins, such as AKT1, MYC, and MAPK1, exhibited significant associations with downregulated DE proteins, including PTEN, TP53, and Caspase3 (Fig. 4A). Similarly, in the 14-DPA samples, upregulated DE proteins, such as TP53, MAPK3, and MYC, displayed notable interactions with top downregulated DE proteins, such as CD4, IL16, and UBC (Fig. 4B).

FIG. 4.

PPI network of top 20 DE proteins potentially regulated by identified circRNAs. PPI network was constructed between upregulated and downregulated proteins in 1-DPA (A) and 14-DPA (B) samples. Darker and brighter gray color show upregulation and downregulation, respectively. PPI, protein–protein interaction.

Data availability statement

All data generated or analyzed during this study are submitted to NCBI SRA databases (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA852678/) and included in supplementary information files.

Discussion

Salamanders are astonishing models due to their remarkable regenerative capacity. Compared with the mammalians' limited renewal potential, axolotl, a well-established salamander species, has been utilized to uncover the molecular and cellular program of regeneration. Despite its evident and significant contribution to extending our understanding of regeneration, the lack of several datasets, such as circRNA data, hampered its complete usage. This study aimed to profile the circRNA composition at the early phase of limb regeneration.

Most of the unique 5331 circRNAs identified in this study were not shared between all time points. Additionally, PCA of circRNA expression highlighted the clustering of samples according to collected time points. These observations may underline circRNA profile alteration during blastema formation phases. Species-specific circRNA among animals (Sun et al., 2019) and tissue-specific circRNA expression within a species (Zhou et al., 2018) were documented before. In line with our data, in a recent study, circRNAs of zebrafish blood, heart, brain, muscle, and gills were compared, and a prominent organ- and tissue-specific circRNA profile was revealed (Sharma et al., 2019). Moreover, a similar origin-specific circRNA expression pattern was detected in human organs (Xu et al., 2017). The profile of circRNAs is more dissimilar for distant organs (Xu et al., 2017) and more alike for the different regions of the same organ, as shown in the brain (Lo et al., 2020).

Although all samples were collected from the limb, the circRNA profile was found considerably distinct among samples, most probably due to the changes in cellular composition and altered biological processes for intact limb (0-DPA), wound healing phase (1-DPA), and blastema stage (14-DPA).

The expression level patterns between examined time points as observed in the RNA-sequencing results are confirmed by qPCR using stably expressed circRNAs with similar expression levels in all time points for normalization. Since there is no established method to assess the function of circRNAs directly, DE circRNAS between 1- and 14-DPA and 0-DPA were queried first to find sponged miRNAs by circRNAs, then the validated target mRNAs of sponged miRNAs.

The regulatory mechanisms of circRNAs on gene expression exhibit a wide range of diversity, and despite the established approaches, the assessment of circRNA's impact on cellular activities lacks a standardized bioinformatics method. The first approach involves linking circRNAs to their respective host genes. Through diverse mechanisms, circRNAs exert regulatory control over their host genes. They function as sponges for RBPs, thereby modulating the stability of these RBPs. Furthermore, circRNAs can interact with RNA polymerase, DNA demethylase, and DNA methyltransferase, influencing the transcriptional regulation of their originating host genes. Consequently, inferring a straightforward regulatory pattern of host genes by circRNAs poses challenges, although several reports have suggested the negative regulation of host genes by circRNAs.

In our dataset, during the initial phase of axolotl limb regeneration, we observed significant downregulation of two circRNAs, namely hsa_circTMEM39B_003 and hsa_circPIP4K2B_011. Notably, a previous study in zebrafish highlighted the involvement of TMEM39B, an endoplasmic reticulum-localized transmembrane protein, in stress resistance and tissue repair (Liu et al., 2022). Additionally, PIP4K2A, a kinase responsible for phosphorylating and removing nuclear PI5P, plays roles in gene expression changes associated with cell-cycle progression, epithelial–mesenchymal transition, reactive oxygen accumulation, and metabolism. A prior investigation underscored its negative regulation of myotube formation (Stijf-Bultsma et al., 2015). Given the crucial roles of TMEM39B in tissue repair and PIP4K2A in suppressing myogenesis, the downregulation of hsa_circTMEM39B_003 and hsa_circPIP4K2B_011 during the early phase of axolotl limb regeneration may augment the expression level and/or stability of TMEM39B and PIP4K2A mRNAs, thereby facilitating a proper regenerative process.

Furthermore, we observed the upregulation of two prominent circRNAs, hsa_circCHKA_005 and hsa_circNCOR1_010. Choline kinase A (CHKA), a key enzyme involved in phosphatidylcholine biosynthesis, is known to be suppressed by a conserved microRNA regulatory circuit during axolotl limb and zebrafish fin regeneration (King and Yin, 2016). Conversely, nuclear receptor corepressor 1 (NCOR1), a chromatin-associated protein, plays a critical role in suppressing fibrosis during axolotl limb regeneration (Rao et al., 2014). Based on existing literature, one would expect the downregulation of CHKA and the upregulation of NCOR1 during axolotl limb regeneration. However, the observed inconsistency may be attributed to the dual regulatory role of circRNAs in both negative and positive regulation of their host genes or the involvement of other gene expression regulatory factors beyond circRNAs.

On the other hand, in the 14-DPA samples, the top-upregulated circRNAs were hsa_circTANC2_010 and hsa_circCHKA_005, while hsa_circEIF4G1_007 and hsa_circPLA2G4C_006 were the top-downregulated circRNAs. TANC2, a multidomain adaptor/scaffolding protein, has a negative regulatory effect on the mTOR pathway, thereby limiting cell survival capacity (Kim et al., 2021). Notably, consistent with the observations in the 1-DPA samples, CHKA was downregulated in the 14-DPA samples. EIF4G, associated with transcriptional elongation, exhibits increased expression levels during axolotl limb regeneration (Voss et al., 2015). Additionally, PLA2G4C has been shown to play a crucial role in tissue repair, particularly in nervous system regeneration (López-Vales et al., 2008). Considering the established roles of these genes and their significance in the regeneration process, as reported in previous studies, we can deduce that the identified circRNAs in this study exert a negative regulatory influence on their respective host genes.

The second approach focuses on the role of circRNAs as miRNA sponges, preventing miRNAs from negatively regulating target mRNAs. Following these approaches, we generated lists of upregulated and downregulated mRNAs and utilized these lists to enrich GO terms, KEGG pathways, and Reactome pathways. In the 1-DPA results, we observed the upregulation of pathways related to telomere length maintenance, cilium assembly, cell cycle regulation, and nitrogen compound transport (Fig. 3A), which aligns well with findings in the literature. Telomere length maintenance is critical for preventing cellular senescence, and it has been documented that axolotls exhibit elevated telomerase activity during limb regeneration to ensure the stability of telomeres for successful renewal (Springhetti et al., 2022). Cilium assembly is initiated in response to developmental signals secreted in the wound zone following axolotl amputation or injury (Wang and Dynlacht, 2018).

Studies have demonstrated that cell cycle-related pathways are activated during limb regeneration, playing a crucial role in blastema formation and regeneration (Sibai et al., 2020; Sibai et al., 2019). Furthermore, a meta-analysis study revealed the upregulation of nitrogen compound transport in blastema tissue during axolotl limb regeneration (Sibai et al., 2019).

Conversely, at 1-DPA, we observed the downregulation of tRNA processing and tRNA modification terms, translation pathways, and Notch signaling-associated pathways. During the early phase of limb regeneration, the downregulation of tRNA processing and modification terms also affects the translation rate. Notch signaling plays a vital role in development and homeostasis. Inhibition of the Notch signaling pathway during the early stages of regeneration has been shown to enhance regeneration capacity in invertebrates and axolotls (El Bejjani and Hammarlund, 2012; Yandulskaya et al., 2022). Hence, the enrichments derived from circRNA analysis are broadly consistent with previously published studies.

In the 14-DPA samples, we observed a remarkable activation of cell cycle-related terms and pathways. During the blastema stage, cells undergo a coordinated process of cell division, which is essential for tissue renewal (Bölük et al., 2022). The accumulation of progenitor cells and their subsequent differentiation into specialized cell types are crucial for proper tissue regeneration, as extensively characterized in previous studies (McCusker et al., 2015; Stocum, 2017). Therefore, the enrichment of cell cycle-related pathways using the circRNA dataset represents a significant finding, highlighting the potential prominent contribution of circRNAs in the regeneration process.

In contrast, those associated with cellular senescence and cell death are noteworthy among the suppressed terms and pathways. Previous studies have demonstrated the elimination of present senescent cells during limb regeneration by macrophages through the axolotl's unique mechanism for sensing and clearing senescent cells (Yun et al., 2015). In line with this, it has been documented that during axolotl limb regeneration, the rate of cell death significantly decreases as a mechanism to contribute to tissue growth following damage or amputation (Wells et al., 2021).

Based on the PPI data, due to the altered levels of circRNAs at 1-DPA, activated key signaling molecules, including, AKT1, MYC, and MAPK1 interacted with suppressed PTEN, TP53, and caspase3. For 14-DPA, there was a significant interaction between upregulated TP53, MAPK3, and MYC, and downregulated CD4, IL16, and UBC. These findings provide significant insights into the regulatory role of circRNAs in the early stages of regeneration process. AKT1, a serine/threonine kinase, plays a crucial role in promoting cell survival, proliferation, and migration, thus facilitating tissue repair (Carbonell et al., 2022). The upregulation of AKT1 suggests its involvement in the activation of pro-regenerative signaling pathways. MYC, a transcription factor, is known for its involvement in cellular growth and proliferation. Its increased expression indicates its contribution to the promotion of cell proliferation during wound healing and blastema phases (Jhamb et al., 2011). MAPK1 is responsible for transmitting extracellular signals to the nucleus and is implicated in cell differentiation and tissue regeneration. Its activation indicates its participation in signaling cascades that promote tissue renewal (Wen et al., 2022).

Conversely, the observed downregulation of PTEN, a tumor suppressor protein, implies its inactivation during wound healing, potentially leading to enhanced cell survival and proliferation (Park et al., 2008).

TP53, a well-known tumor suppressor gene, plays a vital role in regulating cell cycle arrest, DNA repair, and apoptosis. Its downregulation signifies a potential inhibition of apoptotic pathways, enabling cell survival and tissue regeneration (Yun et al., 2013). Additionally, the inactivation of caspase3, a key effector in the apoptotic pathway, further supports the notion of suppressed apoptosis during regeneration (Johnson et al., 2018). CD4, a surface glycoprotein found on T helper cells, plays a critical role in immune system regulation and the activation of immune responses. The observed inactivation of CD4 represents a potential suppression of immune activation during wound healing, which may be important for maintaining an appropriate balance between proinflammatory and anti-inflammatory processes (Li et al., 2020). IL16, an interleukin involved in immune cell recruitment and modulation, is known to contribute to immune responses and inflammation. The inactivation of IL16 indicates a potential dampening of inflammatory responses during wound healing, which could be beneficial for tissue repair and resolution of the inflammatory phase (Mescher, 2017).

Collectively, these findings highlight the intricate regulatory role of circRNAs in modulating the expression and activity of critical signaling molecules involved in wound healing, providing a deeper understanding of the underlying molecular mechanisms in this complex biological process.

The study has several limitations. Human orthologs were used for the annotation and enrichment of the identified circRNAs due to the lack of comprehensive axolotl circRNA and miRNA databases. Additionally, this study lacks wet-lab validations of the identified circRNAs by modulating their expression level during regeneration.

Conclusions

To the best of our knowledge, this study represents the first comprehensive map of circRNAs in axolotl to date. In this study, we present the landscape of circRNAs in axolotl during early phase of limb regeneration. Our investigation describes a list of axolotl circRNAs and reveals that a subset of these circRNAs share orthologous genes with humans, highlighting the potential for translating the messages to human by studying the biological functions and regulatory networks of circRNAs in axolotl. Subsequent bioinformatics analysis indicated that DE circRNAs may contribute to the telomere length maintenance and positively regulate the cell cycle-related processes. Concordantly, circRNAs may suppress cell death and cellular senescence-associated pathways during this process. We anticipate that our findings can serve as a foundation for understanding the circRNAs and their regulatory networks in regeneration, using the axolotl model. These new observations also inform future regenerative medicine research and innovation that harness this model organism.

Footnotes

Authors' Contributions

T.D.: methodology, formal analysis, data curation, visualization, conceptualization, resources, and writing–original draft preparation. B.E.S.: methodology, software, and writing–review and editing. All authors have read and agreed to the final published version of the article.

Author Disclosure Statement

The authors declare they have no conflicting financial interests.

Funding Information

This study was supported by the Scientific and Technological Research Council of Turkey (TÜBİTAK; Project No. 120Z217), and by the BAGEP Award of the Science Academy.

Supplementary Material

Abbreviations Used

References

Ashwal-Fluss

, Meyer

, Pamudurti

, et al. CircRNA biogenesis competes with pre-MRNA splicing. Mol Cell, 2014; 56(1):55–66; doi: 10.1016/j.molcel.2014.08.019

Bölük

, Yavuz

, Demircan

. Axolotl: A resourceful vertebrate model for regeneration and beyond. Dev Dyn, 2022; 251(12):1914–1933; doi: 10.1002/dvdy.520

Camacho

, Coulouris

, Avagyan

, et al. BLAST+: Architecture and applications. BMC Bioinformatics, 2009; 10(1):421; doi: 10.1186/1471-2105-10-421

Carbonell

, Zapata Cardona

, Delgado

. Hydrogen peroxide is necessary during tail regeneration in juvenile axolotl. Dev Dyn, 2022; 251(6):1054–1076; doi: 10.1002/dvdy.386

Demircan

, Hacıbektaşoğlu

, Sibai

, et al. Preclinical molecular signatures of spinal cord functional restoration: Optimizing the metamorphic axolotl (Ambystoma mexicanum) model in regenerative medicine. OMICS J Integr Biol, 2020; 24(6):370–378; doi: 10.1089/omi.2020.0024

Demircan

, İlhan

, Ovezmyradov

, et al. Longitudinal 16S RRNA data derived from limb regenerative tissue samples of axolotl Ambystoma mexicanum. Sci Data, 2019; 6(1):70; doi: 10.1038/s41597-019-0077-7

Demircan

, Keskin

, Dumlu

, et al. Detailed tail proteomic analysis of axolotl (Ambystoma mexicanum) using an mRNA-Seq reference database. Proteomics, 2017; 17(1–2):1600338; doi: 10.1002/pmic.201600338

, Yang

, Liu

, et al. Foxo3 circular RNA retards cell cycle progression via forming ternary complexes with P21 and CDK2. Nucleic Acids Res, 2016; 44(6):2846–2858; doi: 10.1093/nar/gkw027

El Bejjani

, Hammarlund

. Notch signaling inhibits axon regeneration. Neuron, 2012; 73(2):268–278; doi: 10.1016/j.neuron.2011.11.017

10.

Gerber

, Murawala

, Knapp

, et al. Single-cell analysis uncovers convergence of cell identities during axolotl limb regeneration. Science, 2018; 362(6413):eaaq0681; doi: 10.1126/science.aaq0681

11.

Grabherr

, Haas

, Yassour

, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol, 2011; 29(7):644–652; doi: 10.1038/nbt.1883

12.

Hansen

, Jensen

, Clausen

, et al. Natural RNA circles function as efficient microRNA sponges. Nature, 2013a;495(7441):384–388; doi: 10.1038/nature11993

13.

Hansen

, Kjems

, Damgaard

. Circular RNA and miR-7 in cancer. Cancer Res, 2013b;73(18):5609–5612; doi: 10.1158/0008-5472.CAN-13-1568

14.

Jhamb

, Rao

, Milner

, et al. Network based transcription factor analysis of regenerating axolotl limbs. BMC Bioinformatics, 2011; 12(1):80; doi: 10.1186/1471-2105-12-80

15.

Johnson

, Bateman

, DiTommaso

, et al. Systemic cell cycle activation is induced following complex tissue injury in axolotl. Dev Biol, 2018; 433(2):461–472; doi: 10.1016/j.ydbio.2017.07.010

16.

Kim

S-G

, Suho

, Kim

, et al. Tanc2-mediated mTOR inhibition balances mTORC1/2 signaling in the developing mouse brain and human neurons. Nat Commun, 2021; 1(12):2695; doi: 10.1038/s41467-021-22908-4

17.

King

, Yin

. A conserved microRNA regulatory circuit is differentially controlled during limb/appendage regeneration. PLoS One, 2016; 11(6):e0157106; doi: 10.1371/journal.pone.0157106

18.

Kulcheski

, Christoff

, Margis

. Circular RNAs are miRNA sponges and can be used as a new class of biomarker. J Biotechnol, 2016; 238:42–51; doi: 10.1016/j.jbiotec.2016.09.011

19.

Leng

, Dawson

, Thomson

, et al. EBSeq: An empirical Bayes hierarchical model for inference in RNA-Seq experiments. Bioinformatics, 2013; 29(8):1035–1043; doi: 10.1093/bioinformatics/btt087

20.

, Dewey

. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 2011; 12(1):323; doi: 10.1186/1471-2105-12-323

21.

, Liang

, Yang

, et al. Specific ablation of CD4+ T-cells promotes heart regeneration in juvenile mice. Theranostics, 2020; 10(18):8018–8035; doi: 10.7150/thno.42943

22.

, Yang

, Chen

L-L

. The biogenesis, functions, and challenges of circular RNAs. Mol Cell, 2018; 71(3):428–442; doi: 10.1016/j.molcel.2018.06.034

23.

Liu

, Wang

, Shen

, et al. Circbank: A comprehensive database for CircRNA with standard nomenclature. RNA Biol, 2019; 16(7):899–905; doi: 10.1080/15476286.2019.1600395

24.

Liu

, Long

, Liu

, et al. Understanding the function and mechanism of zebrafish Tmem39b in regulating cold resistance. Int J Mol Sci, 2022; 23(19):11442; doi: 10.3390/ijms231911442

25.

Livak

, Schmittgen

. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods, 2001; 25(4):402–408; doi: 10.1006/meth.2001.1262

26.

, Hill

, Vilhjálmsson

, et al. Linking the association between circRNAs and Alzheimer's disease progression by multi-tissue circular RNA characterization. RNA Biol, 2020; 17(12):1789–1797; doi: 10.1080/15476286.2020.1783487

27.

López-Vales

, Navarro

, Shimizu

, et al. Intracellular phospholipase A2 group IVA and group VIA play important roles in Wallerian degeneration and axon regeneration after peripheral nerve injury. Brain, 2008; 131(10):2620–2631; doi: 10.1093/brain/awn188

28.

McCusker

, Bryant

, Gardiner

. The axolotl limb blastema: Cellular and molecular mechanisms driving blastema formation and limb regeneration in tetrapods. Regeneration, 2015; 2(2):54–71; doi: 10.1002/reg2.32

29.

Memczak

, Jens

, Elefsinioti

, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature, 2013; 495(7441):333–338; doi: 10.1038/nature11928

30.

Mescher

. Macrophages and fibroblasts during inflammation and tissue repair in models of organ regeneration. Regeneration, 2017; 4(2):39–53; doi: 10.1002/reg2.77

31.

Pamudurti

, Patop

, Krishnamoorthy

, et al. CircMbl functions in Cis and in trans to regulate gene expression and physiology in a tissue-specific fashion. Cell Rep, 2022; 39(4):110740; doi: 10.1016/j.celrep.2022.110740

32.

Panda

, De

, Grammatikakis

, et al. High-purity circular RNA isolation method (RPAD) reveals vast collection of intronic circRNAs. Nucleic Acids Res, 2017; 45(12):e116; doi: 10.1093/nar/gkx297

33.

Park

, Liu

, Hu

, et al. Promoting axon regeneration in the adult CNS by modulation of the PTEN/MTOR pathway. Science, 2008; 322(5903):963–966; doi: 10.1126/science.1161566

34.

Qin

, Fan

, Wang

, et al. Single-cell RNA-Seq reveals novel mitochondria-related musculoskeletal cell populations during adult axolotl limb regeneration process. Cell Death Differ, 2021; 28(3):1110–1125; doi: 10.1038/s41418-020-00640-8

35.

Rao

, Song

, Jhamb

, et al. Proteomic analysis of fibroblastema formation in regenerating hind limbs of Xenopus laevis froglets and comparison to axolotl. BMC Dev Biol, 2014; 14(1):32; doi: 10.1186/1471-213X-14-32

36.

, Kechris

, Tabakoff

, et al. The MultiMiR R Package and Database: Integration of microRNA–target interactions along with their disease and drug associations. Nucleic Acids Res, 2014; 42(17):e133; doi: 10.1093/nar/gku631

37.

Sharma

, Sehgal

, Mathew

, et al. A genome-wide map of circular RNAs in adult zebrafish. Sci Rep, 2019; 9(1):3432; doi: 10.1038/s41598-019-39977-7

38.

Sibai

, Altuntaş

, Süzek

, et al. Comparison of protein expression profile of limb regeneration between neotenic and metamorphic axolotl. Biochem Biophys Res Commun, 2020; 522(2):428–434; doi: 10.1016/j.bbrc.2019.11.118

39.

Sibai

, Parlayan

, Tuğlu

, et al. Integrative analysis of axolotl gene expression data from regenerative and wound healing limb tissues. Sci Rep, 2019; 9(1):20280; doi: 10.1038/s41598-019-56829-6

40.

Springhetti

, Bucan

, Liebsch

, et al. An identification and characterization of the axolotl (Ambystoma mexicanum, Amex) telomerase reverse transcriptase (Amex TERT). Genes, 2022; 13(2):373; doi: 10.3390/genes13020373

41.

Stijf-Bultsma

, Sommer

, Tauber

, et al. The basal transcription complex component TAF3 transduces changes in nuclear phosphoinositides into transcriptional output. Mol Cell, 2015; 58(3):453–467; doi: 10.1016/j.molcel.2015.03.009

42.

Stocum

. Mechanisms of urodele limb regeneration. Regeneration, 2017; 4(4):159–200; doi: 10.1002/reg2.92

43.

Sun

L-F

, Zhang

, Chen

X-J

, et al. Circular RNAs in human and vertebrate neural retinas. RNA Biol, 2019; 16(6):821–829; doi: 10.1080/15476286.2019.1591034

44.

Szklarczyk

, Gable

, Lyon

, et al. STRING V11: Protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res, 2019; 47(D1):D607–D613; doi: 10.1093/nar/gky1131

45.

Voss

, Palumbo

, Nagarajan

, et al. Gene expression during the first 28 days of axolotl limb regeneration I: Experimental design and global analysis of gene expression. Regeneration, 2015; 2(3):120–136; doi: 10.1002/reg2.37

46.

Vromman

, Yigit

, Verniers

, et al. Validation of circular RNAs using RT-QPCR after effective removal of linear RNAs by ribonuclease R. Curr Protoc, 2021; 1(7):e181; doi: 10.1002/cpz1.181

47.

Wang

, Dynlacht

. The regulation of cilium assembly and disassembly in development and disease. Development, 2018; 145(18):dev151407; doi: 10.1242/dev.151407

48.

Wei

, Fu

, Li

, et al. Single-cell stereo-Seq reveals induced progenitor cells involved in axolotl brain regeneration. Science, 2022; 377(6610):eabp9444; doi: 10.1126/science.abp9444

49.

Wells

, Kelley

, Baumel

, et al. Neural control of growth and size in the axolotl limb regenerate. Elife, 2021; 10:e68584; doi: 10.7554/eLife.68584

50.

Wen

, Jiao

, Tan

. MAPK/ERK pathway as a central regulator in vertebrate organ regeneration. Int J Mol Sci, 2022; 23(3):1464; doi: 10.3390/ijms23031464

51.

Weng

, Wei

, Toden

, et al. Circular RNA CiRS-7—A promising prognostic biomarker and a potential therapeutic target in colorectal cancer. Clin Cancer Res, 2017; 23(14):3918–3928; doi: 10.1158/1078-0432.CCR-16-2541

52.

, Ji

, Zhao

. CircAtlas: An integrated resource of one million highly accurate circular RNAs from 1070 vertebrate transcriptomes. Genome Biol, 2020; 21(1):101; doi: 10.1186/s13059-020-02018-y

53.

Xiao

M-S

, Wilusz

. An improved method for circular RNA purification using RNase R that efficiently removes linear RNAs containing G-Quadruplexes or structured 3′ ends. Nucleic Acids Res, 2019; 47(16):8755–8769; doi: 10.1093/nar/gkz576

54.

, Wu

, Han

, et al. Circular RNA expression profiles and features in human tissues: A study using RNA-Seq data. BMC Genomics, 2017; 18(6):680; doi: 10.1186/s12864-017-4029-3

55.

Yandulskaya

, Miller

, Ansaripour

, et al. Regenerating axolotl retinas regrow diverse cell types with modulation by notch signaling and reconnect to the brain. bioRxiv, 2022; 2022:489898; doi: 10.1101/2022.04.28.489898

56.

, He

Q-Y

. ReactomePA: An R/Bioconductor Package for reactome pathway analysis and visualization. Mol Biosyst, 2016; 12(2):477–479; doi: 10.1039/C5MB00663E

57.

, Wang

L-G

, Han

, et al. ClusterProfiler: An R Package for comparing biological themes among gene clusters. OMICS J Integr Biol, 2012; 16(5):284–287; doi: 10.1089/omi.2011.0118

58.

Yun

, Davaapil

, Brockes

. Recurrent turnover of senescent cells during regeneration of a complex structure. Elife, 2015; 4:e05505; doi: 10.7554/eLife.05505

59.

Yun

, Gates

, Brockes

. Regulation of P53 is critical for vertebrate limb regeneration. Proc Natl Acad Sci U S A, 2013; 110(43):17392–17397; doi: 10.1073/pnas.1310519110

60.

Zhou

, Xie

, Li

, et al. Rat BodyMap transcriptomes reveal unique circular RNA features across tissue types and developmental stages. RNA, 2018; 24(11):1443–1456; doi: 10.1261/rna.067132.118

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.01 MB

2.57 MB

0.01 MB

1.87 MB

0.36 MB

0.06 MB

0.10 MB