Molecular basis of urostyle development in frogs: genes and gene regulation underlying an evolutionary novelty

Abstract

Evolutionary novelties entail the origin of morphologies that enable new functions. These features can arise through changes to gene function and regulation. One key novelty is the fused rod at the end of the vertebral column in anurans, the urostyle. This feature is composed of a coccyx and a hypochord, both of which ossify during metamorphosis. To elucidate the genetic basis of these features, we used laser capture microdissection of these tissues and did RNA-seq and ATAC-seq at three developmental stages in tadpoles of Xenopus tropicalis. RNA-seq reveals that the coccyx and hypochord have two different molecular signatures. Neuronal (TUBB3) and muscle markers (MYH3) are upregulated in coccygeal tissues, whereas T-box genes (TBXT, TBXT.2), corticosteroid stress hormones (CRCH.1) and matrix metallopeptidases (MMP1, MMP8 and MMP13) are upregulated in the hypochord. ATAC-seq reveals potential regulatory regions that are observed in proximity to candidate genes that regulate ossification identified from RNA-seq. Even though an ossifying hypochord is only present in anurans, this ossification between the vertebral column and the notochord resembles a congenital vertebral anomaly seen prenatally in humans caused by an ectopic expression of the TBXT/TBXT.2 gene. This work opens the way to functional studies that can elucidate anuran bauplan evolution.

Headnote

Evolutionary novelties entail the origin of morphologies that enable new functions. These features can arise through changes to gene function and regulation. One key novelty is the fused rod at the end of the vertebral column in anurans, the urostyle. This feature is composed of a coccyx and a hypochord, both of which ossify during metamorphosis. To elucidate the genetic basis of these features, we used laser capture microdissection of these tissues and did RNA-seq and ATAC-seq at three developmental stages in tadpoles of Xenopus tropicalis. RNA-seq reveals that the coccyx and hypochord have two different molecular signatures. Neuronal (TUBB3) and muscle markers (MYH3) are upregulated in coccygeal tissues, whereas T-box genes (TBXT, TBXT.2), corticosteroid stress hormones (CRCH.1) and matrix metallopeptidases (MMP1, MMP8 and MMP13) are upregulated in the hypochord. ATAC-seq reveals potential regulatory regions that are observed in proximity to candidate genes that regulate ossification identified from RNA-seq. Even though an ossifying hypochord is only present in anurans, this ossification between the vertebral column and the notochord resembles a congenital vertebral anomaly seen prenatally in humans caused by an ectopic expression of the TBXT/TBXT.2 gene. This work opens the way to functional studies that can elucidate anuran bauplan evolution.

Keywords:

hypochord, RNA-seq, ATAC-seq, T-box genes, metamorphosis

1. Introduction

Phenotypic and genotypic changes from an ancestral condition undergird the evolution of 'key innovations' [1]. Phenotypic changes of a novel structure reflect changes in the corresponding genotypic/gene regulatory networks [25]. The anuran (frog and toad) urostyle, composed of a coccyx and a hypochord, is morphologically unique from the rest of the vertebrates because of the contribution of an ossifying hypochord, and is therefore considered a structural novelty [6–10]. The coccyx is derived from the paraxial mesoderm, gives rise to the caudal vertebrae [7,11],[12] which subsequently undergo endochondral ossification and fuse together during metamorphosis [6]. The amphibian hypochord, thought to be derived from either endoderm [13–15] or superficial mesoderm [16], is a thin embryonic rod, which degenerates in the rest of anamniotes during early embryonic development but is retained only in frogs and undergoes endochondral ossification during the metamorphic climax [6–10].

Embryonic hypochord in anamniotes is known to have a function in remodelling the dorsal aorta. Previous work [6] has shown how the ossifying hypochord may also play a role in modifying the dorsal aorta by occluding it at the posterior-most end of the hypochord and remodelling it to form two branches, which enter the forelimbs and hind limbs, respectively. Senevirathne et al. [6] speculated that the ossifying hypochord has a role in the evolution of the unique anuran bauplan, which has been conserved in anurans for more than 200 million years [2,17].

The morphological changes of the urostyle are well studied [6,8,9]; however, the molecular mechanisms underlying this unique structure have remained obscure to date. Here, we investigate transcriptomic and gene regulatory networks in the developing urostyle by combining RNA-seq and ATAC-seq approaches. Using previously published morphology work [6] as a framework to identify targeted cells, we used laser capture microdissection (LCM) to reveal the transcriptomics and epigenomics of the two tissue types: coccyx and hypochord.

Paraxial mesoderm-derived coccygeal cells are undifferentiated mesenchymal cells; they undergo chondrification and ossification prior to the initiation of the metamorphic climax [6,7] and could be following a similar gene regulatory network as connective tissues and bones in vertebrates. However, the ossifying hypochord initiates ossification at the onset of metamorphosis. The origin of amphibian hypochordal cells has been hypothesized to be from the endoderm [13–15] or the superficial mesoderm [16]. Regardless of which germ layer it is derived from, hypochord also undergoes endochondral ossification only in anurans, and the hypochord and the coccyx appear to have a conserved regulatory network. Here, we compare the gene expression patterns of coccygeal and hypochordal cells to identify similar/different pathways between the two tissue types, which are derived from two different cell populations.

Using Xenopus tropicalis as our model organism, we address the following questions through this work. Why does the hypochord only ossify in anurans? What are the similarities/differences between the hypochordal and coccygeal molecular pathways? Which genes switch on/off during metamorphosis? By identifying the underlying changes in the genes and gene regulatory networks, our work begins to shed light on the potential genotypic changes underlying a structural novelty.

2. Material and methods

Different stages of X. tropicalis tadpoles were purchased from the National Xenopus Resource (NXR) at the Marine Biological Laboratory (MBL), Woods Hole, MA. Comparisons were made across three significant life-history stages to analyse genes and gene regulation during metamorphosis. The developmental stages used for the experiments were as follows: before metamorphosis/prometamorphic stages (stage 56/57), at the beginning of the metamorphic climax (stage 60/61) and end of metamorphosis (stage 65/66). The tadpoles were euthanized using 0.2% aqueous tricaine methanesulfonate (MS-222), and the specimens were preserved in different fixatives or fresh tissues were taken according to each experiment. Tadpoles were staged according to Nieuwkoop and Faber (NF) [18]. The codes generated for the bioinformatics analyses are deposited in GitHub (https://github.com/GayaniSenevirathne/Senevirathne_et_al_RNAseq.git). The raw sequences are available at NCBI (GSM7701532−GSM7701554) and uploaded to Dryad (https://doi.org/10.5061/dryad.6m905qg68) [19]. Table S1 in the electronic supplementary material contains all the sample names, tissues and replicate details.

Refer to the electronic supplementary material for detailed methodology of RNA and ATAC-seq bioinformatic analyses and HCR in situ hybridization.

2.1. RNA-seq using spatial transcriptomics and laser capture microdissection

All forceps, scissors, surgical blades and lab benches were cleaned with RNAse away and 100% ethanol prior to RNA sequencing experiments. The region where the urostyle forms (demarcated by the 10th and 14th myotomes [6]) was dissected under a Leica L2 light microscope on ice-cold 1× DEPC-treated PBS; all the dissections were done on ice to prevent RNA degradation. The dissected tissue was immediately transferred to ice-cold optimal cutting temperature (OCT) compound and flash frozen in liquid nitrogen and stored at −80°C. For better RNA quality, the tissue blocks were processed the subsequent day.

The two targeted tissue types, coccyx and hypochord, from three individuals at each developmental stage (prometamorphosis, beginning of metamorphic climax and end of metamorphosis) were dissected from frozen sections (figure 1a). Prometamorphic (stage 56) sections of coccyx had undifferentiated mesenchymal cells around the spinal cord, and the hypochord had embryonic hypochordal cells ventral to the notochord. The RNA-sequencing protocol followed a spatial transcriptomics method (Geo-seq; [20]). Prior to sectioning, the cryostat, brushes, adjacent benches/tabletops, blades and pencils/pens were cleaned using RNAse away and 100% ethanol. The frozen tissue blocks were sectioned using a Leica cryostat. The tissue blocks were left inside the cryostat for 20 min, allowing them to equilibrate at −20oC (not doing this resulted in flaky sections or sections breaking when transferred onto the slides). The tissues were sectioned at 10 μM thickness onto PEN membrane 1.0 slides. Five sections were placed on each slide and were allowed to dry at room temperature for 1 min before storing them at −80°C for further processing (the yield of RNA was high when the slides were sectioned on the same day).

On the day of the LCM, slides were thawed at room temperature for 2 min and placed under an UV lamp for 2 min to adhere sections. Next, slides were stained using Cresyl Violet to help visualize cells. For this, slides were graded through an ethanol series. Each wash was 30 s (100% ethanol, 70% ethanol, Cresyl Violet in 70% ethanol and were dehydrated in 70%, 90% and 100% ethanol). Slides were allowed to dry completely before moving to the next step (this step was important to avoid humidity affecting the RNA quality [21]).

The dehydrated slides were processed via LCM with the following settings: aperture (10), speed (20) and energy (50). The hypochordal and coccygeal cells were identified using the ×10 eye piece, and the dissections were done using the ×20. Targeted cells were captured to an adhesive cap Eppendorf tube, with the cap consisting of 50 μl of the lysis buffer. Once cells from coccyx and hypochord were collected (approx. 10 000 cells from 10 sections for each tissue type, 4 replicates were done for each stage and a total of 24 samples), 150 μl of the lysis buffer was added to each tube and the sample was left on ice for 20 min. RNA was extracted from the captured cells using the Takara NucleoSpin RNA XS (cat. no. 740902) kit with slight modification; the filtration step was skipped. cDNA was generated using the SMART-Seq v4 Ultra Low Input RNA Kit for Sequencing with the number of amplification cycles set to 18. cDNA was purified using Agencourt Ampure XP magnetic beads (Beckman Coulter) and was sequenced using the HiSeq PE100. Downstream bioinformatic RNA-seq analysis is given in detail under electronic supplementary material.

2.2. Gene regulation and ATAC-seq

The same developmental stages that were used for the RNA-seq studies were analysed. The urostyle region was dissected out as a fresh chunk of tissue. The OMNI-ATAC-seq protocol was used to identify open chromatin regions during urostyle development. Two replicates from each developmental stage, coinciding with the RNA-seq and morphological studies, were selected.

The tadpoles were anaesthetized in MS-222, dissected on ice-cold 1× PBS, and mechanically crushed using a pestle that had been cleaned using 100% ethanol prior to this step to obtain a homogenized sample. These steps were done on ice to prevent protein degradation. After homogenization, cells were counted using the BioRad Tc20 automated cell counter. All samples consisted of 75 000−100 000 cells. The subsequent steps followed the OMNI ATAC-seq protocol [22,23] with slight modifications using the Illumina Tagment DNA Enzyme and Buffer kit: cells were lysed in an ice-cold lysis buffer, followed by a transposition step using Tn5 Transposase, and DNA was purified using the Zymo DNA Clean and Concentrator. Purified DNA was amplified with 13 amplification cycles (the number of cycles was optimized by an additional qPCR step). Finally, the libraries were purified using the Zymo DNA Clean and Concentrator and were sequenced using a NovaSeq 2000 (100 bp paired end). The ATAC-seq bioinformatic analysis is discussed in detail under electronic supplementary material.

3. Results and discussion

3.1. Disparity in gene expression profiles of the coccyx and hypochord

The hypochord and coccyx show considerable differences in cellular composition and differentiation (figures 1a and 2), and the gene expression profiles directly reflect this (figures 1b,c). The two tissue types fuse at the end of metamorphosis, coinciding with the degeneration of the notochord. The total analysis consisted of 21 458 genes; both tissue types and the three timepoints were used as factors in the DESeq2 analysis where a binomial generalized linear model was implemented. Principal component analysis (PCA) revealed that the coccygeal and hypochordal samples generate two separate clusters (figure 1b), and a heatmap showed the two tissue types possess two different gene expression profiles (figure 1c). A total of 3298 genes were differentially expressed between the urostyle and hypochord across development (the FDR < 0.05), among these DEGs, 2828 genes were significantly upregulated and 470 were downregulated in the coccygeal region compared hypochord.

During coccyx development, the spinal cord and peripheral nerves remodel, and muscular organization changes. Consistent with these anatomical changes, we find that there are DEGs involved in differentiation and development of the nervous system (e.g. NEUROD6, PRDM12, COCH and APBA2 [24,25]), genes that are expressed during skeletal muscle development (e.g. ACTN2 [26]) and genes that are directly involved in chondrocyte/osteocyte differentiation (e.g. RUNX2, COL9A1 and SOX8) (figure 3) [27,28].

3.2. Transcriptomic differences across different time points during urostyle development

There are numerous studies of metamorphic transcriptomes (e.g. [29–35]) but none on the urostyle. We first looked into urostyle-responsive transcriptomes by comparing genes that are differentially expressed in the coccyx and hypochord at different time points: Before metamorphosis versus beginning of metamorphosis (electronic supplementary material, S1B–S3) and before/beginning of metamorphosis versus end of metamorphosis (electronic supplementary material, figure S1A). This analysis identified 5664 DEGs that fell within the thresholds of FDR < 0.01 (adjusted p-values < 0.05 and log fold change of 1.5) and showed unique expression patterns that were significant at each time point.

Several sets of genes were upregulated and downregulated across the three developmental time points (figure 1c). We identified four unique clusters when the transcriptomes were compared between the three developmental time points (before and beginning of metamorphosis versus end of metamorphosis (electronic supplementary material, figure S1). Cluster A has 47 genes that were highly downregulated at the end of metamorphosis (switched off) compared with the other two time points. This cluster includes genes involved in muscle contraction and M-band stabilization in fast skeletal muscles (e.g. TRDN and MYOM2I [36,37]), skeletal development (e.g. SOX9 [38]), response to inflammation (PTX3 [39]), filament organizing genes (e.g. KRT18.I and VIM.2 [40,41]), extracellular matrix organizing and connective tissue-strengthening (e.g. COL9A1, COL8A1 and CHAD [42,43]) and stress regulation (CRCH.1 [44]). The other two gene clusters, B and C (electronic supplementary material, figure S1A), comprise genes that are both downregulated and upregulated at the end of metamorphosis. Cluster C also has 15 genes that are downregulated at the end of metamorphosis, which include collagen markers (e.g. COL9A3) and skeletal muscle function genes (e.g. MYL1 and ACTN3 [45,46]). Genes that are upregulated (10 genes) are within Cluster B and are involved in mitosis (CCNB1 [47]), development of neurons (POU3F1 [48]) and maintenance of myelin sheath (PLP1 [49]). When before metamorphosis was compared with beginning/end of metamorphosis, clustering of the 100 top-most significant genes revealed metamorphic genes that were downregulated before metamorphosis but were upregulated during metamorphosis. Heatmap clustering revealed five main clusters (electronic supplementary material, figure S1B). Cluster A included 28 genes that were downregulated (switched off) before metamorphosis in both coccyx and hypochord, but as soon as metamorphosis was initiated, these genes were upregulated; they are involved in functions like collagen synthesis (SERPINH1 [50]), cell cycle (CDK6 [51]) and thyroid hormone inactivation (DIO3; [52]). Clusters B and C include genes that are switched on prior to metamorphosis and are switched off at the onset of metamorphosis: HES8, FOXP2, EGR1, HOXD11 and PVALB are representative examples. Cluster D is enriched with genes that are involved in blood sugar control (e.g. THRAP3, IGF2BP3 [53,54]), which are downregulated before metamorphosis but are upregulated at the onset of metamorphosis. This part of the transcriptomic analysis identified DEGs that are specific to the three significant time points (before metamorphosis versus onset of metamorphosis versus end of metamorphic climax).

We next explored the Gene Ontology (GO) function of these significant genes during development. The DEGs and the corresponding p-values from the differential expression analyses were imported into an online database of reactome pathways ('Reactome pathway browser') to compare the functional aspect of these genes (figure 3a; electronic supplementary material,

figure S2). DEGs upregulated before metamorphosis were enriched for GO terms like 'DNA replication and pre-initiation',

'synthesis of DNA', 'polymerase switching', 'G1/S transition' (figure 3a,b), whereas the DEGs upregulated during metamorphosis include genes that function in 'collagen formation', 'cross linking of collagen fibrils', 'RUNX2 regulated bone development' and 'osteocyte differentiation' (figures 3a,c,d).

3.3. Hypochord, metamorphosis and T-box genes

As there are no data on the genes that are expressed during hypochordal ossification, we used the DEGs identified by the coccyx versus hypochord comparisons to scrutinize this. This identified 470 genes that were uniquely upregulated only within the hypochordal tissues (they fell within the significant threshold of adjusted p value < 0.05 and FDR < 0.01). We saw that hypochordal tissues express high concentrations of TBX1, TBXT, TBXT.2 and HAND2 (figure 3c).

T-box genes have been implicated in early mesodermal patterning and, especially, Brachyury/Xbra/TBXT is essential in early mesodermal formation [55–59], and TBXT homologues across vertebrates induce the mesoderm [60,61]. T-box genes have already been identified as being pivotal components in the differentiation of the posterior axial column (e.g. [55–58,62–65]) and seem to be playing a role in hypochordal ossification as well. Xenopus has two paralogues of the gene Brachyury: TBXT (also known as Xbra or T) and TBXT.2 (also known as Xbra3 or T2). When Brachyury is knocked out, it causes loss of posterior mesoderm and failure to differentiate the notochord [63,66]. However, the expression of TBXT and TBXT.2 in late-developing tadpole structures has not been reported so far.

The temporal and spatial expression patterns of TBXT and TBXT.2 make them good candidate genes for regulating ossification only in hypochordal tissues. We performed HCR in situ hybridization to examine the temporal and spatial expression patterns of TBXT in hypochordal ossification. TBXT expression is exclusively concentrated along the ossifying hypochord at the onset of metamorphosis but is not evident in prometamorphic nor at the end of metamorphic climatic tadpoles (electronic supplementary material, figure S6).

The pelvic region undergoes dramatic changes during metamorphosis, and this period is thought to represent the developmental stage that is most susceptible to predation. The underlying stress of the remodelling tissues and hormonal responses can also be seen by the increased expression of CRCH.1 (corticosol steroid stress hormones), having a normal hormonal response to stress. Other than these genes, the hypochord also expresses significant concentrations of VEGF and HAND2. These two genes are involved in vascular development and can also be seen expressed in embryonic hypochord where VEGF plays a role in the formation of the hypochord (e.g. [13,14,67]).

3.4. Transcriptomic comparisons between coccyx + hypochord and other ossifying elements

Ossification happens through two major processes: endochondral (cartilaginous precursors used as a template) and intramembranous (direct ossification of the condensed mesenchymal cells) [68]. Even though coccyx and hypochord are highlighted as derivatives of two different cell populations, they both undergo endochondral ossification [6]. During this process, mesenchymal cells condense (commit to form osteoprogenitors) and aggregate to form cartilaginous precursors during early development. Cartilaginous precursors expand and cells proliferate, next the extracellular matrix is synthesized, and finally, mineralization of the matrix occurs. These steps are similar to other bones in vertebrates, which undergo endochondral ossification as well [69].

Youlten et al. [70] identified three clusters of GO functions during osteocyte development: (i) an 'early expression cluster' (expressed in osteoprogenitors/osteoblast-like cells); (ii) an 'early activation cluster' (expressed in early osteocytes); and (iii) a 'maturation cluster' (expressed in mature osteocytes). We compared the spatial and temporal transcriptomic maps of the osteocytes introduced by Youlten et al. [70] with the coccygeal and hypochordal transcriptomics to see whether the molecular underpinning of ossification is similar in the genes responsible for the formation of the urostyle as well.

Early expression cluster. This included GO term functions 'extracellular matrix organization', 'angiogenesis', 'cartilage development' and 'connective tissue development' (electronic supplementary material, figure S4). Out of the genes that are differentially expressed, there are some that are inactive before metamorphosis in the hypochord (e.g. COL22A1, COL16A1, COL6A3, RUNX1 and IHH) but are highly expressed once the metamorphosis is initiated. High expression of these genes in the coccygeal cells even before the onset of metamorphosis corroborates our morphological studies, where we revealed that the post caudal vertebrae of the coccyx initiated mesenchymal cell aggregation early in development (1.5 months after embryogenesis) versus 2 months in hypochord. Apart from the differences in the temporal expression of genes within the 'early expression cluster', a few genes involved in cartilage development are not present in the hypochord compared with the coccyx (e.g. FOXL1, RUNX3, FOXD3, PMM2 and EDN1).

Early activation cluster. This cluster includes the GO terms 'axon guidance', 'axon development', 'axogenesis', 'regulation of axogenesis' and 'neuron projection guidance' (figure 4). While the coccyx DEGs act in a similar way to the rest of the long bones in vertebrates within this cluster, hypochord shows a different pattern. Most of the genes (e.g. NTRN, SLITRK3, POUF42 and DCC) that are discussed as essential regulators in guiding the axons in long bones are not expressed within the hypochord (figure 4).

Maturation cluster. The GO term functions 'bone development', 'skeletal system development', 'regulation of ossification', 'ossification' and 'osteoblast differentiation' are included in this cluster (electronic supplementary material, figure S5). Maturation period in the hypochord happens once the metamorphosis is initiated and when the tadpole reaches the end of its metamorphic climax (electronic supplementary material, figure S5). Within the hypochord, genes involved in ossification (e.g. GPC3, TMEM19, IFITM5, COL11A1, PHOSPHO1 and SOX8) and osteoblast differentiation (e.g. GLI1, FBN2 and SATB2) are highly expressed in tadpoles at the end of the metamorphic climatic and are inactive at prometamorphic stages. Comparatively, in the coccyx, since the ossification happens prior to the metamorphic climax, the majority of the genes are highly expressed even at the beginning of metamorphosis. A few genes (e.g. TBX15, BARX2, SHH and AXIN2) are not expressed in hypochord nor in the coccyx, compared with the other ossifying long bones in vertebrates.

This transcriptomic comparison led to three main findings (i) Between the two tissue types, the coccyx's DEGs share similarities with the other bones' transcriptomics in vertebrates. (ii) Hypochord undergoes its early activation period before metamorphosis, and a maturation period once metamorphosis is initiated. (iii) Hypochordal DEGs lack an early activation period, which includes most of the axon developing genes.

3.5. ATAC-seq and urostyle-responsive gene regulation

During anuran metamorphosis, the larval body form undergoes dramatic remodelling within 6–8 days, and this is reflected in both morphological and gene expression patterns. Therefore, it can be extrapolated that gene regulation changes over this same time period. To study the underlying changes in chromatin accessibility, we used an ATAC-seq approach (figure 5a) using the same developmental stages and the same number of replicates as the RNA-seq work. The number of peaks varied between the three stages that we used (electronic supplementary material). More than 50% of peaks were distributed in distal intergenic regions. The rest of the peaks were distributed along intronic, exons and promoter regions. When comparing the three time points, the most significant change of peak distribution observed was the percentage of peaks that fell on the exon regions (other than the first exon): before metamorphosis, the percentage was lower (<1%) when compared with the number of peaks that were seen at the beginning and at the end of metamorphosis (7−10%) (figure 5b).

Next, we compared the ATAC-seq data with the RNA-seq data and observed that majority of the peaks are located close to the upregulated genes in the hypochord and coccyx that were identified from the transcriptomic data. The gene TBXT, which is upregulated in the hypochord, has peaks located within the exon 1, 2 and 3 before metamorphosis, and the exon 1's peak is lost during metamorphosis (figure 5c). Other genes expressed in hypochordal tissues like MMP1 also show a similar pattern where the open chromatin region changes during the course of metamorphosis (figure 5c). Genes that were upregulated in the coccyx, e.g. TUBB3 shows significant peaks that are present throughout development (figure 5c). This could be because the coccygeal ossification occurs early in development (after 1.5 months) compared with the hypochord, which is initiated at the onset of metamorphosis. These results highlight urostyle-responsive regulatory regions during development and need further scrutinization using functional assays.

4. Conclusion

The coccyx and hypochord both undergo endochondral ossification during development, and similar ossification patterns were reflected in the gene expression profiles as well. Though there were major differences in some transcriptomes (e.g. presence of T-box genes, CRCH.1 and MMPs in hypochordal tissues at the onset of metamorphosis versus absent in the coccyx), there were similarities in genes that were involved in endochondral ossification: we show that genes that are involved in cartilage and bone formation, extracellular matrix organization and thyroid hormone responsive elements are present in both tissues (electronic supplementary material, figure S2) but differ temporally (the coccyx starts ossifying after 1.5 months, whereas the hypochord initiates its ossification only at the onset of metamorphosis).

4.1. T-box genes and the hypochord

The coccyx and hypochord have two sets of differentially expressed genes. This analysis revealed a large set of genes (electronic supplementary material) that are uniquely upregulated in the hypochord and have not been reported before in the context of a developing hypochord. One of the most significant groups of genes that is upregulated in the hypochord is the T-box genes (TBXT and TBXT.2). T-box genes have a 180 bp DNA-binding domain that is highly conserved. Orthologues of the gene Brachyury, one of the highly expressed T-box genes in the hypochord, are present in all multicellular organisms [62]. Brachyury is important in posterior mesoderm development (initially expressed in the developing mesoderm but later restricted to the tail bud and notochord) [65]. While early mesoderm differentiation patterning depends highly on TBXT/TBXT.2, a role for these genes in later developmental stages has not been previously reported or discussed. We hypothesize that the presence of high levels of TBXT/TBXT.2 causes the undifferentiated hypochordal cells to undergo ossification at the onset of metamorphosis. Such unusual ossification appears to also occur in response to a congenital vertebral column malformation (VCM) in humans that happens because of a Brachyury gene mutation in the intron 7 [64] and in the highly conserved T-box sequence [71]; these VCMs eventually lead to sacral agenesis (frog-like) syndrome in babies. In humans with this abnormality, increased expression or duplications of the TBXT gene result in production of excess Brachyury [62,72]. Apart from these mutations, TBXT/TBXT.2 genes also induce EMT in humans when over expressed in carcinoma cells [73], and it has also been recorded that duplications of the Brachyury gene cause vertebral column chordomas [73].

Previous studies have shown that Brachyury acts as a switch in posterior mesoderm specification during embryogenesis and is restricted to the anteroposterior axis [55,56]. Here, during hypochordal ossification, the onset of metamorphosis could be triggering ectopic expression of TBXT/TBXT.2 in hypochordal cells, which could potentially express posterior mesodermal genes and subsequently activate downstream targets of TBXT/TBXT2, which in turn initiates chondrification and ossification. When Brachyury genes are highly expressed in human chordoma cells, matrix metalloproteinases (e.g. MMP12, MMP13 and MMP24; 74) are also upregulated at the same time (which is also seen in hypochordal cells). The extent to which the human and frog conditions are similar awaits functional tests.

4.2. Coccyx and hypochord versus other vertebrate skeletal elements

The coccyx and hypochord undergo endochondral ossification and show an array of genes that are similar to the genes expressed in other long bones that undergo endochondral ossification in vertebrates (e.g. mesenchymal-to-chondrocytes involved genes like BMPs, SOX9; chondrocytes-to-osteoblasts/osteocytes were seen in highly expressed genes like RUNX2, Osterix and IHH). Apart from these similarities, when comparing the already published osteocyte transcriptomics [70], hypochord shows some considerable differences among the rest of the bones in vertebrates. Hypochordal cells express osteoprogenitor-specific genes before the metamorphic climax, and metamorphosis acts as a switch that activates osteogenesis (versus in coccyx osteogenesis is initiated prior to metamorphosis). Other than the temporal differences observed regarding ossification, the DEGs of the hypochord reveal that hypochordal cells lack the 'early activation phase', which includes regulators needed in 'Axogenesis' and 'Axon development' in ossifying bones (figure 4). Vertebrate bones are innervated by sensory and sympathetic nerves during skeleton development [75], where the periosteum and bone marrow have the highest density of nerves whereas the mineralized matrix has very few [75–77]. During development, bone innervation and endochondral ossification happen simultaneously [75], and it is hypothesized that axon guidance regulates formation of the neuronal network, which is subsequently required for the osteocyte network formation [70]. It is surprising that the ossifying hypochord lacks the genes needed for axon development (figure 4), and our results raise the possibility that the hypochordal development may be disconnected from the neuronal signals. Future work is needed scrutinizing the innervation patterns within the hypochord during its development to better understand this.

Our integrative approach, using morphological and molecular datasets (genes and gene regulation) on the development of the urostyle, scrutinizes the evolution of a novelty. Future work targeting the candidate genes responsible for the development of the urostyle, together with functional assays, will shed light on the evolution of this structural enigma.

Ethics. All experiments were approved by the Marine Biological Laboratory (Protocol 18-42) and the University of Chicago (Protocol 72615).

Data accessibility. The raw sequences are available at NCBI (GSM7701532-GSM7701554) and uploaded to Dryad [78]. Table S1 of the electronic supplementary material contains all the sample names, tissues and replicate details [79].

Declaration of AI use. We have not used AI-assisted technologies in creating this article.

Authors' contributions. G.S.: conceptualization, formal analysis, methodology, validation, visualization, writing-original draft, writing-review and editing; N.H.S.: conceptualization, resources, supervision, writing-review and editing.

Both authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration. We declare we have no competing interests.

Funding. This work was supported by the University of Chicago Biological Sciences and the Brinson Foundation (to N.H.S.) and by the O'Brien and Hasten Fellowship to G.S.

Acknowledgements. We would like to thank Marko Horb, Nikko-Ideen Shaidaini and Marcin Wlizla of the National Xenopus Resource (NXR), Marine Biological Laboratory, for husbandry and providing X. tropicalis tadpoles. We also thank Thomas Stewart, James Hanken, Victoria Prince, Madhava Meegaskumbura and Shubin Lab members for helpful discussions on this work.

Sidebar

Cite this article: Senevirathne G, Shubin NH. 2024 Molecular basis of urostyle development in frogs: genes and gene regulation underlying an evolutionary novelty. Open Biol. 14: 240111.

Received: 30 April 2024

Accepted: 25 July 2024

Subject Areas:

genetics, developmental biology, molecular biology

Author for correspondence:

Gayani Senevirathne

e-mail: [email protected]

References

1. Galis F, Metz JAJ. 2007 Evolutionary novelties: the making and breaking of pleiotropic constraints. Integr. Comp. Biol. 47, 409–419. (doi:10.1093/icb/icm081)

2. Shubin N, Tabin C, Carroll S. 2009 Deep homology and the origins of evolutionary novelty. Nature 457, 818–823. (doi:10.1038/nature07891)

3. Tarazona OA, Slota LA, Lopez DH, Zhang G, Cohn MJ. 2016 The genetic program for cartilage development has deep homology within bilateria. Nature 533, 86–89. (doi:10.1038/ nature17398)

4. Tschopp P, Tabin CJ. 2017 Deep homology in the age of next-generation sequencing. Phil. Trans. R. Soc. B 372, 20150475. (doi:10.1098/rstb.2015.0475)

5. Wagner GP. 2015 Evolutionary innovations and novelties: let us get down to business! Zool. Anz. 256, 75–81. (doi:10.1016/j.jcz.2015.04.006)

6. Senevirathne G, Baumgart S, Shubin N, Hanken J, Shubin NH. 2020 Ontogeny of the anuran urostyle and the developmental context of evolutionary novelty. Proc. Natl Acad. Sci. USA 117, 3034–3044. (doi:10.1073/pnas.1917506117)

7. Handrigan GR, Wassersug RJ. 2007 The anuran bauplan: a review of the adaptive, developmental, and genetic underpinnings of frog and tadpole morphology. Biol. Rev. Camb. Philos. Soc. 82, 1–25. (doi:10.1111/j.1469-185X.2006.00001.x)

8. Branham AE, List JC. 1979 Development of the urostyle during metamorphosis in five species of anurans. J. Morphol. 159, 311–329. (doi:10.1002/jmor.1051590303)

9. Kovalenko EE, Danilov IG. 2006 Diversity of the sacral-urostyle region in the family Bufonidae (Amphibia, Anura). 1. Actual diversity of sacrum in Bufonidae. Zool. Zh. 85, 500–516.

10. Snell CA. 2015 Identifying ranid urostyle, ilial and anomolous bones from a 15th century London well. Herpetol. J. 25, 245–255.

11. Sánchez RS, Sánchez SS. 2013 Characterization of pax1, pax9, and uncx sclerotomal genes during Xenopus laevis embryogenesis. Dev. Dyn. 242, 572–579. (doi:10.1002/dvdy.23945)

12. Sánchez RS, Sánchez SS. 2015 Paraxis is required for somite morphogenesis and differentiation in Xenopus laevis. Dev. Dyn. 244, 973–987. (doi:10.1002/dvdy.24294)

13. Cleaver O, Krieg PA. 1998 VEGF mediates angioblast migration during development of the dorsal aorta in Xenopus. Development 125, 3905–3914. (doi:10.1242/dev.125.19.3905)

14. Cleaver O, Seufert DW, Krieg PA. 2000 Endoderm patterning by the notochord: development of the hypochord in Xenopus. Development 127, 869–879. (doi:10.1242/dev.127.4.869)

15. Löfberg J, Collazo A. 1997 Hypochord, an enigmatic embryonic structure: study of the axolotl embryo. J. Morphol. 232, 57–66. (doi:10.1002/(SICI)1097-4687(199704)232:1<57: AID-JMOR3>3.0.CO;2-L)

16. Shook DR, Majer C, Keller R. 2004 Pattern and morphogenesis of presumptive superficial mesoderm in two closely related species, Xenopus laevis and Xenopus tropicalis. Dev. Biol.270, 163–185. (doi:10.1016/j.ydbio.2004.02.021)

17. Shubin NH, Jenkins FA. 1995 An early jurassic jumping frog. Nature 377, 49–52. (doi:10.1038/377049a0)

18. Nieuwkoop PD, Faber J. 1994 Normal table of. X. laevis 252, 598.

19. Senevirathne G, Shubin N. Molecular basis of urostyle development: genes and gene regulation underlying an evolutionary novelty. Dryad Digital Repository (doi:10.5061/dryad. 6m905qg68)

20. Chen J, Suo SB, Tam PPL, Han JDJ, Peng GD, Jing NH. 2017 Spatial transcriptomic analysis of cryosectioned tissue samples with geo-seq. Nat. Protoc. 12, 566–580. (doi:10.1038/ nprot.2017.003)

21. Ordway GA, Szebeni A, Duffourc MM, Dessus-Babus S, Szebeni K. 2009 Gene expression analyses of neurons, astrocytes, and oligodendrocytes isolated by laser capture microdissection from human brain: detrimental effects of laboratory humidity. J. Neurosci. Res. 87, 2430–2438. (doi:10.1002/jnr.22078)

22. Buenrostro JD, Wu B, Chang HY, Greenleaf WJ. 2015 ATAC-seq: a method for assaying chromatin accessibility genome-wide. Curr. Protoc. Mol. Biol. 109, 21. (doi:10.1002/ 0471142727.mb2129s109)

23. Corces MR et al. 2017 An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues. Nat. Methods 14, 959–962. (doi:10.1038/nmeth.4396)

24. Uittenbogaard M, Baxter KK, Chiaramello A. 2010 NeuroD6 genomic signature bridging neuronal differentiation to survival via the molecular chaperone network. J. Neurosci. Res. 88, 33–54. (doi:10.1002/jnr.22182)

25. Rahman MM, Kim IS, Ahn D, Tae HJ, Park BY. 2020 PR domaincontaining protein 12 (prdm12) is a downstream target of the transcription factor zic1 during cellular differentiation in the central nervous system: PR domain containing protein is the right form. Int. J. Dev. Neurosci. 80, 528–537. (doi:10.1002/jdn.10048)

26. Mills M, Yang N, Weinberger R, Vander Woude DL, Beggs AH, Easteal S, North K. 2001 Differential expression of the actin-binding proteins, alpha-actinin-2 and -3, in different species: implications for the evolution of functional redundancy. Hum. Mol. Genet. 10, 1335–1346. (doi:10.1093/hmg/10.13.1335)

27. Youlten SE et al. 2020 Osteocyte transcriptome mapping identifies a molecular landscape controlling skeletal homeostasis and susceptibility to skeletal disease. bioRxiv. (doi:10. 1101/2020.04.20.051409)

28. Qin X et al. 2020 Runx2 is essential for the transdifferentiation of chondrocytes into osteoblasts. PLoS Genet. 16, e1009169. (doi:10.1371/journal.pgen.1009169)

29. Brown DD, Cai L. 2007 Amphibian metamorphosis. Dev. Biol.306, 20–33. (doi:10.1016/j.ydbio.2007.03.021)

30. Callery EM, Elinson RP. 2000 Thyroid hormone-dependent metamorphosis in a direct developing frog. Proc. Natl Acad. Sci. USA 97, 2615–2620. (doi:10.1073/pnas.050501097)

31. Zhao L, Liu L, Wang S, Wang H, Jiang J. 2016 Transcriptome profiles of metamorphosis in the ornamented pygmy frog Microhyla fissipes clarify the functions of thyroid hormone receptors in metamorphosis. Sci. Rep. 6, 27310. (doi:10.1038/srep27310)

32. Wang S, Liu L, Liu J, Zhu W, Tanizaki Y, Fu L, Bao L, Shi YB, Jiang J. 2019 Gene expression program underlying tail resorption during thyroid hormone-dependent metamorphosis of the ornamented pygmy frog Microhyla fissipesFront. Endocrinol. (Lausanne) 10, 11. (doi:10.3389/fendo.2019.00011)

33. Yaoita Y, Brown DD. 1990 A correlation of thyroid hormone receptor gene expression with amphibian metamorphosis. Genes Dev. 4, 1917–1924. (doi:10.1101/gad.4.11.1917)

34. Kanamori A, Brown DD. 1996 The analysis of complex developmental programmes: amphibian metamorphosis. Genes Cells 1, 429–435. (doi:10.1046/j.1365-2443.1996.d01-251.x)

35. Brown DD, Wang Z, Kanamori A, Eliceiri B, Furlow JD, Schwartzman R. 1995 Amphibian metamorphosis: a complex program of gene expression changes controlled by the thyroid hormone. Recent Prog. Horm. Res. 50, 309–315. (doi:10.1016/b978-0-12-571150-0.50018-4)

36. Giacomazzi G et al. 2017 MicroRNAs promote skeletal muscle differentiation of mesodermal ipsc-derived progenitors. Nat. Commun. 8, 1249. (doi:10.1038/s41467-017-01359-w)

37. Auxerre-Plantié E et al. 2020 Identification of MYOM2 as a candidate gene in hypertrophic cardiomyopathy and tetralogy of fallot, and its functional evaluation in the drosophila heart. Dis. Model. Mech. 13, dmm045377. (doi:10.1242/dmm.045377)

38. Hattori T et al. 2010 SOX9 is a major negative regulator of cartilage vascularization, bone marrow formation and endochondral ossification. Development 137, 901–911. (doi:10. 1242/dev.045203)

39. Magrini E, Mantovani A, Garlanda C. 2016 The dual complexity of PTX3 in health and disease: a balancing act? Trends Mol. Med. 22, 497–510. (doi:10.1016/j.molmed.2016.04.007)

40. Velez-delValle C, Marsch-Moreno M, Castro-Muñozledo F, Galván-Mendoza IJ, Kuri-Harcuch W. 2016 Epithelial cell migration requires the interaction between the vimentin and keratin intermediate filaments. Sci. Rep. 6, 24389. (doi:10.1038/srep24389)

41. Gan Z et al. 2016 Vimentin intermediate filaments template microtubule networks to enhance persistence in cell polarity and directed migration. Cell Syst. 3, 252–263.(doi:10.1016/ j.cels.2016.08.007)

42. Brachvogel B et al. 2013 Comparative proteomic analysis of normal and collagen IX null mouse cartilage reveals altered extracellular matrix composition and novel components of the collagen IX interactome. J. Biol. Chem. 288, 13481–13492. (doi:10.1074/jbc.M112.444810)

43. Hessle L et al. 2014 The skeletal phenotype of chondroadherin deficient mice. PLoS One 8, e63080. (doi:10.1371/journal.pone.0063080)

44. Reul JMHM, Holsboer F. 2002 On the role of corticotropin-releasing hormone receptors in anxiety and depression. Dialogues Clin. Neurosci. 4, 31–46. (doi:10.31887/DCNS.2002.4.1/ jreul)

45. Schiaffino S, Rossi AC, Smerdu V, Leinwand LA, Reggiani C. 2015 Developmental myosins: expression patterns and functional significance. Skelet. Muscle 5, 22. (doi:10.1186/s13395015-0046-6)

46. Pickering C, Kiely J. 2017 ACTN3: more than just a gene for speed. Front. Physiol. 8, 1080. (doi:10.3389/fphys.2017.01080)

47. Strauss B, Harrison A, Coelho PA, Yata K, Zernicka-Goetz M, Pines J. 2018 Cyclin B1 is essential for mitosis in mouse embryos, and its nuclear export sets the time for mitosis. J. Cell Biol. 217, 179–193. (doi:10.1083/jcb.201612147)

48. Zhu Q et al. 2014 The transcription factor Pou3f1 promotes neural fate commitment via activation of neural lineage genes and inhibition of external signaling pathways. elife 3, e02224. (doi:10.7554/eLife.02224)

49. Gould RM, Oakley T, Goldstone JV, Dugas JC, Brady ST, Gow A. 2008 Myelin sheaths are formed with proteins that originated in vertebrate lineages. Neuron Glia Biol. 4, 137–152. (doi:10.1017/S1740925X09990238)

50. Widmer C, Gebauer JM, Brunstein E, Rosenbaum S, Zaucke F, Drögemüller C, Leeb T, Baumann U. 2012 Molecular basis for the action of the collagen-specific chaperone hsp47/ SERPINH1 and its structure-specific client recognition. Proc. Natl Acad. Sci. USA 109, 13243–13247. (doi:10.1073/pnas.1208072109)

51. Tigan AS, Bellutti F, Kollmann K, Tebb G, Sexl V. 2016 CDK6-a review of the past and a glimpse into the future: from cell-cycle control to transcriptional regulation. Oncogene 35, 3083–3091. (doi:10.1038/onc.2015.407)

52. Bianco AC, da Conceição RR. 2018 The deiodinase trio and thyroid hormone signaling. Methods Mol. Biol. 1801, 67–83. (doi:10.1007/978-1-4939-7902-8_8)

53. Choi JH et al. 2014 Thrap3 docks on phosphoserine 273 of PPAR gamma and controls diabetic gene programming. Genes Dev. 28, 2361–2369. (doi:10.1101/gad.249367.114)

54. Dong H, Zhang Y, Wang J, Kim DS, Wu H, Sjögren B, Gao W, Luttrell L, Wang H. 2017 Regulator of G protein signaling 2 is a key regulator of pancreatic β-cell mass and function. Cell Death Dis. 8, e2821. (doi:10.1038/cddis.2016.216)

55. Cunliffe V, Smith JC. 1994 Specification of mesodermal pattern in Xenopus laevis by interactions between brachyury, noggin and xwnt-8. EMBO J. 13, 349–359. (doi:10.1002/j. 1460-2075.1994.tb06268.x)

56. Cunliffe V, Smith JC. 1992 Ectopic mesoderm formation in Xenopus embryos caused by widespread expression of a brachyury homologue. Nature 358, 427–430. (doi:10.1038/ 358427a0)

57. Hayata T, Eisaki A, Kuroda H, Asashima M. 1999 Expression of brachyury-like T-box transcription factor, Xbra3 in Xenopus embryo. Dev. Genes Evol. 209, 560–563. (doi:10.1007/ s004270050289)

58. Messenger NJ, Kabitschke C, Andrews R, Grimmer D, Núñez Miguel R, Blundell TL, Smith JC, Wardle FC. 2005 Functional specificity of the Xenopus T-domain protein brachyury is conferred by its ability to interact with smad1. Dev. Cell 8, 599–610. (doi:10.1016/j.devcel.2005.03.001)

59. Smith JC, Price BM, Green JB, Weigel D, Herrmann BG. 1991 Expression of a xenopus homolog of brachyury (T) is an immediateearly response to mesoderm induction. Cell 67, 7987. (doi:10.1016/0092-8674(91)90573-h)

60. Schulte-Merker S, Smith JC. 1995 Mesoderm formation in response to brachyury requires FGF signalling. Curr. Biol. 5, 62–67.

61. Yasuoka Y, Shinzato C, Satoh N. 2016 The mesoderm-forming gene brachyury regulates ectoderm-endoderm demarcation in the coral acropora digitifera. Curr. Biol. 26, 2885–2892. (doi:10.1016/j.cub.2016.08.011)

62. Chen M, Wu Y, Zhang H, Li S, Zhou J, Shen J. 2020 The roles of embryonic transcription factor BRACHYURY in tumorigenesis and progression. Front. Oncol. 10, 961. (doi:10.3389/ fonc.2020.00961)

63. Gentsch GE, Spruce T, Monteiro RS, Owens NDL, Martin SR, Smith JC. 2018 Innate immune response and off-target mis-splicing are common morpholino-induced side effects in Xenopus. Dev. Cell 44, 597–610.(doi:10.1016/j.devcel.2018.01.022)

64. Ghebranious N et al. 2008 A missense T (brachyury) mutation contributes to vertebral malformations. J. Bone Miner. Res. 23, 1576–1583. (doi:10.1359/jbmr.080503)

65. Hotta K, Takahashi H, Asakura T, Saitoh B, Takatori N, Satou Y, Satoh N. 2000 Characterization of brachyury-downstream notochord genes in the Ciona intestinalis embryo. Dev. Biol.224, 69–80. (doi:10.1006/dbio.2000.9765)

66. Paraiso KD, Blitz IL, Zhou JJ, Cho KWY. 2019 Morpholinos do not elicit an innate immune response during early xenopus embryogenesis. Dev. Cell 49, 643–650.(doi:10.1016/j.devcel. 2019.04.019)

67. Cleaver O, Tonissen KF, Saha MS, Krieg PA. 1997 Neovascularization of the Xenopus embryo. Dev. Dyn. 210, 66–77. (doi:10.1002/(SICI)1097-0177(199709)210:1<66:AID-AJA7>3. 0.CO;2-#)

68. Breeland G, Sinkler MA, Menezes RG. 2021 Embryology, bone ossification. Treasure Island, FL: StatPearls.

69. Mackie EJ, Ahmed YA, Tatarczuch L, Chen KS, Mirams M. 2008 Endochondral ossification: how cartilage is converted into bone in the developing skeleton. Int. J. Biochem. Cell Biol. 40, 46–62. (doi:10.1016/j.biocel.2007.06.009)

70. Youlten SE et al. 2021 Osteocyte transcriptome mapping identifies a molecular landscape controlling skeletal homeostasis and susceptibility to skeletal disease. Nat. Commun. 12, 2444. (doi:10.1038/s41467-021-22517-1)

71. Postma AV et al. 2014 Mutations in the T (brachyury) gene cause a novel syndrome consisting of sacral agenesis, abnormal ossification of the vertebral bodies and a persistent notochordal canal. J. Med. Genet. 51, 90–97. (doi:10.1136/jmedgenet-2013-102001)

72. Zhu J, Kwan KM, Mackem S. 2016 Putative oncogene brachyury (T) is essential to specify cell fate but dispensable for notochord progenitor proliferation and EMT. Proc. Natl Acad. Sci. USA 113, 3820–3825. (doi:10.1073/pnas.1601252113)

73. Henderson SR et al. 2005 A molecular map of mesenchymal tumors. Genome Biol. 6, R76. (doi:10.1186/gb-2005-6-9-r76)

74. Wan Z, Jiang D, Chen S, Jiao J, Ji L, Shah AS et al. 2016 T-box transcription factor brachyury promotes tumor cell invasion and metastasis in non-small cell lung cancer via upregulation of matrix metalloproteinase 12. Oncol. Rep. 36, 306–314.

75. Tomlinson RE, Christiansen BA, Giannone AA, Genetos DC. 2020 The role of nerves in skeletal development, adaptation, and aging. Front. Endocrinol. (Lausanne) 11, 646. (doi:10. 3389/fendo.2020.00646)

76. Mach DB et al. 2002 Origins of skeletal pain: sensory and sympathetic innervation of the mouse femur. Neuroscience 113, 155–166. (doi:10.1016/s0306-4522(02)00165-3)

77. Castañeda-Corral G, Jimenez-Andrade JM, Bloom AP, Taylor RN, Mantyh WG, Kaczmarska MJ, Ghilardi JR, Mantyh PW. 2011 The majority of myelinated and unmyelinated sensory nerve fibers that innervate bone express the tropomyosin receptor kinase A. Neuroscience 178, 196–207. (doi:10.1016/j.neuroscience.2011.01.039)

78. Koshkina V et al. 2018 Data from: Integrated species distribution models: combining presence-background data and site-occupany data with imperfect detection. Dryad Digital Repository. (doi:10.5061/dryad.8467g)

79. Senevirathne G, Shubin N. 2024Data from: Molecular basis of urostyle development in frogs: genes and gene regulation underlying an evolutionary novelty. Figshare. (doi:10.6084/ m9.figshare.c.7410576)