Temporal transcriptomic profiling reveals dynamic changes in gene expression of Xenopus animal cap upon activin treatment

Activin, a member of the transforming growth factor-β (TGF-β) superfamily of proteins, induces various tissues from the amphibian presumptive ectoderm, called animal cap explants (ACs) in vitro. However, it remains unclear how and to what extent the resulting cells recapitulate in vivo development. To comprehensively understand whether the molecular dynamics during activin-induced ACs differentiation reflect the normal development, we performed time-course transcriptome profiling of Xenopus ACs treated with 50 ng/mL of activin A, which predominantly induced dorsal mesoderm. The number of differentially expressed genes (DEGs) in response to activin A increased over time, and totally 9857 upregulated and 6663 downregulated DEGs were detected. 1861 common upregulated DEGs among all Post_activin samples included several Spemann’s organizer genes. In addition, the temporal transcriptomes were clearly classified into four distinct groups in correspondence with specific features, reflecting stepwise differentiation into mesoderm derivatives, and a decline in the regulation of nuclear envelop and golgi. From the set of early responsive genes, we also identified the suppressor of cytokine signaling 3 (socs3) as a novel activin A-inducible gene. Our transcriptome data provide a framework to elucidate the transcriptional dynamics of activin-driven AC differentiation, reflecting the molecular characteristics of early normal embryogenesis.


Results
Validation of mesoderm-inducing activity by activin A treatment. We first validated the effective concentrations of activin A for the induction of dorsal mesoderm genes by semiquantitative RT-PCR. As shown in Supplementary Fig. S1, 50 ng/mL activin A predominantly induced dorsal mesoderm genes, compared with 10 ng/mL activin A, which is consistent with previous findings that higher concentrations of activin induce dorsal mesoderm [11][12][13] . We then validated the ability of activin A at 50 ng/mL concentration to induce dorsal mesoderm at later time points. As shown in Supplementary Fig. S2A, activin A-treated ACs exhibited extensive elongation, which is reportedly caused by the differentiation of notochord and somite 10 . Furthermore, we found that the expression of the notochord marker, chrd, and the somite markers, myf5, actc1 and myogenin were induced after 48 h of activin A treatment ( Supplementary Fig. S2B). These data indicate that 50 ng/mL activin A is the optimal concentration for ACs differentiation into dorsal mesoderm.

Transcriptomic changes in activin A-treated ACs.
To comprehensively understand the temporal gene expression during activin A-induced ACs differentiation, we collected ACs immediately after dissection from blastula (Pre_activin; sibling embryos reached stages 8. [5][6][7][8][9], and ACs after the cultivation in activin A solution for 1, 3, 6 and 9 h (Post 1h_activin; stage 9.5, Post 3h_activin; stages 10-10.25, Post 6h_activin; stages 10.5-11, and Post 9h_activin; stage 11.5, respectively) (Fig. 1A), and performed a time-course RNA-seq to analyze the transcriptome. To assess the dynamics of gene expression variability following activin A treatment, we performed principal component analysis (PCA) and hierarchical clustering analysis. PCA using all samples showed that there was a clear distinction between the same groups at five time points, and the transcriptomic changes in ACs from Pre_activin to Post 9h_activin (Fig. 1B). Moreover, in hierarchical clustering analysis, we found that Pre_activin samples were distinguished from Post_activin groups (Post 1 h, 3 h, 6 h, and 9h_activin), which were further divided into two subgroups, short-term (Post 1 h and 3h_activin) and long-term activin A treatment (Post 6 h and 9h_activin) (Fig. 1C). These data suggest that the gene expression of ACs is significantly changed after activin A treatment in a temporal manner.

Differentially expressed genes profiles between Pre_activin and Post_activin samples.
We next analyzed the statistically significant differentially expressed genes (DEGs) between activin A-treated and -untreated samples. Based on a false discovery rate (FDR) of less than 0.05 and a log2 fold change (FC) ≥ 2 or ≤ − 2, we extracted DEGs in Post_activin samples compared to the Pre_activin sample as a control. The full list of DEGs is shown in Supplementary Table S1, and the top five upregulated and downregulated DEGs in each comparison are shown in Table 1. As shown in Fig. 2A, volcano plots showed that 2573, 4149, 6323, and 8858 genes were upregulated, and 111, 1901, 4704 and 6008 genes were downregulated in the Post 1 h, 3 h, 6 h, and 9h_activin samples, respectively. Thus, the number of DEGs gradually increased over time after treatment with activin A, but more upregulated DEGs were detected than downregulated DEGs across all time points (Fig. 2B). As shown in Fig. 2C, the ratio of unique DEGs (non-overlapping DEGs among different Post_activin samples) in Post 9h_activin sample was higher than other Post_activin samples, implying that time point-or context-specific gene expression actively occurs at least 9 h after treatment with activin A.
According to the data set of DEGs shown in the Venn diagrams, 1861 genes including activin-responsive targets, mix1, gsc, eomes, foxa4 and nodal1 [32][33][34][35][36] , were commonly upregulated ( Fig. 3A  Classification of genes that exhibited distinct expression pattern. To identify the set of genes with similar temporal expression patterns, we performed k-means clustering for 2000 genes, which were variably expressed after activin A treatment. We identified four clusters with distinct expression patterns ( Fig. 4A and Supplementary Table S3). Cluster A (967 genes) contained a set of genes that was downregulated over time after activin A treatment. Clusters B, C, and D contained genes that were upregulated after activin A treatment, but there was little difference in gene expression patterns. Cluster B (213 genes) contained genes whose expression was significantly upregulated at 1 or 3 h and subsequently downregulated at later time points. Cluster C (204 genes) contained genes whose expression was markedly from 3 or 6 h, and slightly downregulated but maintained at higher expression levels at later time points. Cluster D (616 genes) contained genes that were drastically upregulated at 6 or 9 h. To compare activin A-induced ACs differentiation in vitro with in vivo development, we presented the temporal gene expression of five representatives of each cluster in activin A-treated ACs and that in whole gastrula embryos. The transcriptome dataset of the whole gastrula embryos was obtained from published RNA-seq data 25 . As shown in Fig. 4B, the downregulation of genes in cluster A was observed in both cultured ACs and whole embryos. Conversely, the expression of genes in cluster B was actively induced in ACs at the early time point (in Post 3 h and 6h_activin ACs), but these genes were modestly induced in whole embryos (Fig. 4B). The expression of genes in cluster C and cluster D increased in ACs later (in Post 9h_activin ACs), and similar temporal changes in gene expression were observed in whole embryos (Fig. 4B).
To further investigate the functional diversity of DEGs over time after activin A treatment, we next performed functional enrichment analysis on each cluster. The top ten GO terms in each cluster are listed in Supplementary Table S4. GO enrichment analysis of biological processes showed that genes in cluster A were associated with nuclear envelop organization ( Table S4). In addition, the predicted protein-protein interaction (PPI) network for genes in cluster C showed that some genes, encoding transcription factors such as OTX2, SOX17, T, and LHX1, were assigned in the mesoderm development GO term (GO: 0007498) ( Fig. 4C pink circles, Supplemental Table S5).

Upregulation of socs3 in response to activin A. The most significantly upregulated DEGs in activin
A-treated ACs at the early time point was socs3 (Table 1, Post 1h_activin). Therefore, we focused on socs3 as a possible activin-inducible gene. socs3 is reportedly involved in cytokine signaling, and induced in response to wound healing after epithelium trauma 37 . To determine whether socs3 is upregulated in response to activin A treatment, we examined the temporal expression of socs3 in ACs with or without activin A treatment using semiquantitative RT-PCR. As shown in Fig. 5, socs3 was slightly expressed in ACs before the cultivation (− activin, 0 h), but significantly upregulated after 1 and 3 h of activin A treatment (+ activin, 1 h and 3 h), compared with in the absence of activin A treatment at the same time points (− activin, 1 h and 3 h). These data suggest that socs3 was maternally expressed or induced in response to the traumatic dissection from blastulae, but further induced after activin A treatment. Thus, we validated that RNA-seq can be used to identify novel activin-inducible genes in ACs. To predict the socs3-associated genes in activin A-treated ACs, we inferred a gene network using tempo- Table 1. List of top five upregulated and downregulated DEGs. DEGs were identified in the comparison of each Post_activin vs. Pre_activin sample. Sequences that did not show the significant hits or annotations (listed as model IDs starting with "loc" and "Xelaev") were represented as undescribed 1-8. X. laevis is an allotetraploid, and contains two sets of genes, L and S genes.

Discussion
In this study, we captured the transcriptional features of ACs that had the potential to develop into dorsal mesoderm. Although the transcriptional organization in activin A-treated ACs did not completely correspond to DMZ-expressing genes, activin A-treated ACs expressed a variety of dorsal mesoderm genes that are involved in Spemann's organizer formation. We also found that temporal gene expression patterns were similar between activin A-treated ACs and whole embryos, with some notable differences. We also provided the data on expression of socs3 as a potential activin-inducible gene, which was identified from our RNA-seq dataset.
Activin signaling regulates other signaling pathways for mesoderm development in Xenopus. For example, activin reportedly induces the expression of fgf3/8, and also induces the negative regulator of FGF signaling, dusp1, for the dorsoventral patterning of mesoderm 40 . Furthermore, activin signaling reportedly activates the expression of downstream genes of p53 signaling, p21waf1 and pal-1, and the mesendoderm gene, mix.2, through   Table S4). Thus, activin signaling can modulate other signaling pathways for establishment of mesoderm in ACs. Notably, the expression of siamois (sia1 and sia2) that acts downstream of Wnt signaling and is essential for the formation of the Spemann's organizer 44 was not significantly induced in activin A-treated ACs (Supplementary Table S7). Furthermore, all previously reported DMZ-enriched genes did not overlap with upregulated DEGs in response to activin A (Fig. 3C). These differences in gene expression between differentiated ACs and embryonic development are possibly due to the limited effects of activin A at this concentration in ACs, or due to other signaling derived from the neighboring cells in the embryo, but are not in ACs. www.nature.com/scientificreports/ In this study, we classified upregulated genes whose expression reached the peak at early (Post 1 h and 3 h), intermediate (Post 3 h and 6 h), and later time points (Post 6 h and 9 h) following activin A treatment ( Fig. 4A; corresponding to cluster B, C and D). The expression of activin-early responsive genes (cluster B) was actively induced in ACs than in whole gastrula embryos (Fig. 4B). In addition, the temporal trends in expression of intermediate-or late-responsive genes (clusters C and D) in ACs were highly consistent with those observed in gastrulae, although there were some differences in transcription levels and temporal association (Fig. 4B). These data imply that the primary or direct effects of activin A on early-responsive genes result in the downstream gene expression in concordance with normal development. As for the activin-responsive transcription, several activin-responsive targets were identified from the experiments using the combinatorial treatment of activin and the protein synthesis inhibitor, cycloheximide 26,27,33 . In addition, genes associated with transcription factors of activin/nodal signaling, Smad2/3 and FoxhI, were determined by the chromatin immunoprecipitation sequencing (ChIP-seq) using X. tropicalis gastrulae 45,46 . We confirmed that several activin-responsive target genes and Smad2/3-FoxhI associated genes were upregulated in ACs at early time points after activin A treatment (Supplementary Table S3 and Supplementary Fig. S8). It suggests the possibility that the effects of activin A on direct target genes lead to the induction of indirect or secondary responsive genes, resulting in diverse gene expressions in ACs later.
Interestingly, numerous genes previously undescribed in mesoderm formation were observed. Especially, downregulated genes in ACs following activin A treatment have not been investigated well in former experiments, as compared to upregulated genes ( Table 1). We first demonstrated that genes associated with the organization of nuclear envelop were downregulated in ACs following activin A treatment, for example, ctdnep1 and lemd3 ( Fig. 4C and Supplementary Table S4). CTDNEP1, encoding a serine/threonine protein phosphatase, functions as a mediator of lipid composition in nuclear envelop, but it represses both TGF-β and BMP signal transduction [47][48][49][50] . Lemd3, encoding a nuclear membrane protein, also suppresses the TGF-β signal transduction through physical interaction with downstream transcription factors of TGF-β signaling, Smads 51,52 . These data suggest the possibility that the nuclear envelop is organized during ACs differentiation into the endodermal and mesodermal lineage, leading to reinforcement of the TGF-β signal activation. The role of downregulated genes associated with the regulation of golgi and endoplasmic reticulum and rRNA/tRNA transcription in early developmental processes remains to be elucidated (Fig. 4C), but our data may lead to the identification of novel molecular mechanisms in vivo development.
During the blastula stage of embryonic development, zygotic gene activation occurs, which is called the midblastula transition (MBT) 53,54 . It has been reported that the open chromatin accessibility for gene expression becomes higher after MBT (at least after stage 9) in Xenopus ACs 55 . Our transcriptome data using ACs dissected from blastula demonstrated the dynamic expression of various epigenetic regulators, which increase chromatin accessibility and lead to acceleration of gene transcription. For example, various hdac genes, encoding histone deacetylases, were downregulated (except for hdac9.S, hdac10.L, and hdac7), whereas ep300 and hat1, encoding histone acetyltransferase, were upregulated in response to activin A treatment (Supplementary Tables S1 and  S2). In general, histone acetylation targeted to lysines progressively increases transcriptional activation, whereas HDAC deacetylation activity represses transcription 56 . Our data demonstrating downregulation of hdac and upregulation of ep300 and hat suggest the possibility that gene expression is activated by epigenetic mechanisms in ACs following activin A treatment. Another example of epigenetic regulators that alter the expression levels in response to activin A treatment is dnmt1, encoding DNA methyltransferase (Supplementary Tables S1 and S2). During normal development, Dnmt1 suppresses premature gene activation before MBT 57 , and the expression levels are downregulated from late blastula to early gastrula for zygotic gene activation 58 . Consistent with the data, we detected the downregulation of dnmt1 in ACs upon activin A treatment (Supplementary Tables S1 and S2). Taken together, these data imply that the in vitro culture system using ACs treated with activin A recapitulates zygotic gene activation controlled by epigenetic mechanisms in normal development.
We also newly found that socs3 was upregulated by activin A treatment (Table 1 and Fig. 5). socs3 plays a role in skin inflammation 59 , and it has been reported to control the balance between stem cell pluripotency and differentiation into the primitive endoderm lineage via FGF-MAPK pathway 60,61 . In the inferred network, a pluripotency factor, pou5f3.2.S, an endoderm marker, gata4.S, and a downstream kinase of MAPK, map3k1.S, were obtained as potential socs3-associated genes (Supplementary Fig. S6). These previous findings with the predicted socs3 subnetwork indicated that socs3 is involved in the maintenance of ACs pluripotency and activin A-driven differentiation. Further research is needed to elucidate whether activin signaling regulates the pluripotent and differentiation states in Xenopus ACs via the modulation of socs3 expression.
Taken together, our results demonstrating the developmental potency of ACs relied on activin A provide information required for early development. Our transcriptome in a mesendoderm-biased environment may serve as a resource for understanding the network during primary germ layer formation in combination with in vivo transcriptome datasets [23][24][25] and previous GRN resources [62][63][64][65] . The in vitro ACs assay has some limitations, such as fewer cells and limiting growth factors at local levels, but it enables to understand features that could not be captured in complex processes within embryos.

Methods
Preparation of Xenopus embryos. The wildtype X. laevis strains were maintained in accordance with guidelines on animal housing 66 . All animal procedures were approved by the Animal Experimentation Committee of the Teikyo University (approval numbers: 18-020), and authors carried out all methods in our experiments in accordance with ARRIVE guidelines. Embryos were artificially fertilized, dejellied, and incubated in 0.1 × Steinberg's solution 67 . Embryos were cultured until the blastula stages (stages 8. [5][6][7][8][9]   Reverse transcription PCR (RT-PCR) assay. Equal amounts of total RNA extracted from each sample were subjected to RT-PCR as described previously 71 . PCR was performed using an EmeraldAmp PCR Master kit (Takara) using the following protocol: 28 cycles of 98 °C for 10 s, 55 °C for 30 s, and 72 °C for 30 s, followed by final elongation at 72 °C for 10 min. The sequences of the tested primers are listed in Supplementary Table S8. Electrophoresis with 1% agarose gel was carried out to verify the band size of the PCR products. For semiquantitative analysis, the band intensity of PCR products was quantified using Bio Image Intelligent Quantifier software (Bio Image Systems, Inc.) and the relative expression of the tested genes was calculated by dividing the band intensity of the tested genes with that of eef1α1 as an internal control. Three pools in each sample were assayed for RT-PCR. Experiments were carried out at least twice, and a representative result is shown when similar results were obtained.

Functional enrichment analysis and the inferred network. Gene ontology (GO) analysis and PPI
network analysis were performed after conversion of the X. laevis genes to the human ortholog genes as previously described 72 . GO analysis and KEGG pathway analysis (Supplementary Table S4 Use of published datasets. The dorsal marginal zone (DMZ)-enriched genes were recovered from Kakebeen et al. 23 with an additional step to indicate the enrichment of DMZ at stage 11: log transformed counts (DMZ)-log transformed counts (VMZ) > 0.5. The DMZ-and ventral marginal zone (VMZ)-enriched genes were referred to Ding et al. 24 . The gene expression datasets in the whole embryos were referred to Session et al. 25 . Smad2/3-associated genes were recovered from stage 10.5 peak datasets in Gupta et al. 45 , and Smad2/3 and FoxhI target genes from Chiu et al. 46 .