==== Front PLoS One PLoS One plos PLOS ONE 1932-6203 Public Library of Science San Francisco, CA USA 10.1371/journal.pone.0287858 PONE-D-23-06337 Research Article Biology and Life Sciences Developmental Biology Metamorphosis Biology and Life Sciences Developmental Biology Life Cycles Tadpoles Biology and Life Sciences Genetics Gene Expression Gene Regulation Biology and Life Sciences Genetics Gene Expression Biology and Life Sciences Cell Biology Cell Processes Cell Cycle and Cell Division Biology and Life Sciences Cell Biology Chromosome Biology Chromatin Biology and Life Sciences Genetics Epigenetics Chromatin Biology and Life Sciences Genetics Gene Expression Chromatin Biology and Life Sciences Genetics Genomics Biology and Life Sciences Biochemistry Hormones Thyroid Hormones Metamorphic gene regulation programs in Xenopus tropicalis tadpole brain Molecular basis for metamorphosis of tadpole brain Raj Samhitha Conceptualization Data curation Formal analysis Investigation Validation Writing – original draft Writing – review & editing 1 Sifuentes Christopher J. Conceptualization Data curation Formal analysis Investigation Software Validation Visualization Writing – original draft Writing – review & editing 1 ¤a https://orcid.org/0000-0001-8472-6332 Kyono Yasuhiro Conceptualization Data curation Formal analysis Investigation Writing – review & editing 2 ¤b https://orcid.org/0000-0002-2773-8174 Denver Robert J. Conceptualization Data curation Formal analysis Funding acquisition Project administration Supervision Visualization Writing – original draft Writing – review & editing 1 2 * 1 Department of Molecular, Cellular and Developmental Biology, University of Michigan, Ann Arbor, Michigan, United States of America 2 Neuroscience Graduate Program, University of Michigan, Ann Arbor, Michigan, United States of America Coen Laurent Editor Museum National d’Histoire Naturelle, FRANCE Competing Interests: The authors have declared that no competing interests exist. ¤a Current address: Chan Zuckerberg Initiative, Redwood City, CA, United States of America ¤b Current address: Department of Pharmacy Practice, College of Pharmacy, University of Illinois at Chicago, Chicago, IL, United States of America * E-mail: rdenver@umich.edu 29 6 2023 2023 18 6 e02878583 3 2023 14 6 2023 © 2023 Raj et al 2023 Raj et al https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Amphibian metamorphosis is controlled by thyroid hormone (TH), which binds TH receptors (TRs) to regulate gene expression programs that underlie morphogenesis. Gene expression screens using tissues from premetamorphic tadpoles treated with TH identified some TH target genes, but few studies have analyzed genome-wide changes in gene regulation during spontaneous metamorphosis. We analyzed RNA sequencing data at four developmental stages from the beginning to the end of spontaneous metamorphosis, conducted on the neuroendocrine centers of Xenopus tropicalis tadpole brain. We also conducted chromatin immunoprecipitation sequencing (ChIP-seq) for TRs, and we compared gene expression changes during metamorphosis with those induced by exogenous TH. The mRNA levels of 26% of protein coding genes changed during metamorphosis; about half were upregulated and half downregulated. Twenty four percent of genes whose mRNA levels changed during metamorphosis had TR ChIP-seq peaks. Genes involved with neural cell differentiation, cell physiology, synaptogenesis and cell-cell signaling were upregulated, while genes involved with cell cycle, protein synthesis, and neural stem/progenitor cell homeostasis were downregulated. There is a shift from building neural structures early in the metamorphic process, to the differentiation and maturation of neural cells and neural signaling pathways characteristic of the adult frog brain. Only half of the genes modulated by treatment of premetamorphic tadpoles with TH for 16 h changed expression during metamorphosis; these represented 33% of the genes whose mRNA levels changed during metamorphosis. Taken together, our results provide a foundation for understanding the molecular basis for metamorphosis of tadpole brain, and they highlight potential caveats for interpreting gene regulation changes in premetamorphic tadpoles induced by exogenous TH. http://dx.doi.org/10.13039/100000154 Division of Integrative Organismal Systems 1456115 https://orcid.org/0000-0002-2773-8174 Denver Robert J. This research was supported by National Science Foundation grant IOS 1456115 to RJD. https://www.nsf.gov/ The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Data AvailabilityThe RNA-seq datasets used for the current analyses were previously published and the data deposited in GEO (developmental: GSE139267; 3,5,3’-triiodothyronine (T3)-induced: GSE130816). The TR ChIP-seq dataset has been deposited in GEO (GSE214697). Data Availability The RNA-seq datasets used for the current analyses were previously published and the data deposited in GEO (developmental: GSE139267; 3,5,3’-triiodothyronine (T3)-induced: GSE130816). The TR ChIP-seq dataset has been deposited in GEO (GSE214697). ==== Body pmcIntroduction Most amphibian species have complex life cycles with two discrete life history stages separated by a metamorphosis. Amphibian metamorphosis is controlled by thyroid hormone (TH), which binds to ligand-activated transcription factors (TH receptors—TRs) to regulate gene transcription [1]. All vertebrates studied to date have two genes that code for TRs, designated alpha and beta. The TRs bind to DNA as heterodimers with retinoid X receptor at TH response elements (TREs), which are comprised of two, six nucleotide stretches called half sites; the most common type of TRE has two half sites in direct repeat with a four base spacer (DR+4) [2]. The TH-TR complexes induce biochemical, physiological and morphological changes during metamorphosis by controlling transcription of genes encoding proteins that underlie tissue transformations [3]. Thyroid hormone biosynthesis and release is regulated by neuroendocrine centers in tadpole brain whose maturation depends on TH [4,5]. Neurosecretory neurons located in the anterior preoptic area (POa) and hypothalamus produce neurohormones that are released into the pituitary portal system where they travel to the adenohypophysis to regulate secretion of thyroid stimulating hormone, which then acts on thyroid follicle cells to induce TH biosynthesis and release [4]. To understand the molecular basis for tadpole metamorphosis, investigators have used gene expression screening techniques like subtractive hybridization, differential display, DNA microarray, and recently, RNA sequencing to identify genes regulated by TH in tadpole tissues. These experiments identified genes that encode proteins involved with TH signaling like TRs, monodeiodinases, chromatin modifiers, and transcription factors that mediate TH actions, among others [6–24]. While these studies have illuminated some of the cellular pathways regulated by TH, and have provided sets of genes for further study, most relied on treatment of premetamorphic tadpoles with TH versus analysis of spontaneous metamorphosis. Therefore, we currently lack a comprehensive understanding of the molecular changes that drive the transformation of tadpole tissues during spontaneous metamorphosis. We recently published an RNA-seq dataset generated from Xenopus tropicalis tadpole brain (region of the preoptic area/thalamus/hypothalamus) at four developmental stages encompassing the beginning to the end of spontaneous metamorphosis [25]. In that study we investigated the relationship between DNA methylation and gene regulation changes during metamorphosis, but we did not analyze cellular pathways represented in the RNA-seq dataset. Ta and colleagues [24] recently reported an RNA sequencing experiment conducted on the midbrain of Xenopus laevis tadpoles at four developmental stages from premetamorphosis to just before the climax of metamorphosis. To our knowledge, there have been no comprehensive, genome-wide analyses of gene expression changes in any tadpole tissue that spanned the start to the end of spontaneous metamorphosis, nor has there been a comparison of these changes with those induced by exogenous TH. Such analyses are essential to understand the gene regulation programs during spontaneous metamorphosis, and to determine if the gene expression changes induced by exogenous TH in premetamorphic tadpoles reflect what occurs during normal development. In the current study, we investigated gene regulation changes in neuroendocrine centers of X. tropicalis tadpole brain during spontaneous metamorphosis. We analyzed RNA-seq data [25] collected at four stages from the beginning to the completion of metamorphosis. We also conducted chromatin precipitation sequencing (ChIP-seq) for TR on tadpole brain chromatin at metamorphic climax when circulating TH and expression of TRs are maximal. Lastly, we compared gene regulation changes during spontaneous metamorphosis to the set of genes modulated by treatment of premetamorphic tadpoles with 3,5,3’-triiodothyronine (T3) for 16 h [21]. Our findings provide essential data for understanding the molecular basis for metamorphosis of tadpole brain, and also demonstrate potential caveats for interpreting data from TH-induced gene regulation studies conducted on premetamorphic tadpoles. Materials and methods Animal care and use We obtained X. tropicalis tadpoles by in-house breeding, reared them in dechlorinated tap water (water temperature 25°C, pH 7) and maintained them on a 13L:11D photoperiod. We fed tadpoles ad libitum with pulverized frog brittle (NASCO, Fort Atkinson, WI) or Sera Micron plankton food. Tadpoles were staged using the developmental staging system of Nieuwkoop and Faber [26] (NF). Animals were euthanized following humane methods to alleviate suffering by immersion in 0.1% benzocaine followed by decapitation. All procedures involving animals were conducted under an approved animal use protocol (PRO00006809) in accordance with the guidelines of the Institutional Animal Care and Use Committee at the University of Michigan. RNA sequencing (RNA-seq) The RNA-seq datasets used for the current analyses were previously published [21,25] and the data deposited in GEO (developmental: GSE139267; 3,5,3’-triiodothyronine (T3)-induced: GSE130816). For the developmental gene expression analysis, we measured mRNA levels in tadpole preoptic area/thalamus/hypothalamus (S1 Fig) at four stages of metamorphosis (NF50: premetamorphosis, when the larvae grow but little or no morphological change occurs and plasma TH concentrations are low; NF56: prometamorphosis, when hindlimb growth accelerates and plasma TH concentration rises; NF62: metamorphic climax, the most rapid phase of morphological change when thyroid activity is at its peak; NF66: completion of metamorphosis, juvenile frog; n = 3/developmental stage) [25]. For analysis of T3-induced gene expression changes we measured mRNA by RNA-seq in preoptic area/thalamus/hypothalamus of premetamorphic tadpoles (NF54) treated with or without T3 (5 nM dissolved in the aquarium water) for 16 h (n = 3/treatment) [21]. We re-analyzed the two RNA-seq datasets with the same bioinformatics tools as before [21,25], but using the most recent X. tropicalis genome build (v9.1) from Ensembl and setting a false discovery rate (FDR) adjusted p value <0.05. We generated heatmaps for the top 100 upregulated and the top 100 downregulated genes determined during the developmental interval NF50 to NF62 using the software Heatmapper [27]. We clustered the genes in the developmental RNA-seq dataset into groups based on their expression profiles using the software Clust [28]. We used gProfiler software [29] to conduct functional enrichment analysis of gene ontology (GO), KEGG pathway, and reactome (REACT) pathways within the set of differentially expressed genes for each comparison. KEGG pathway data were rendered using Pathview software [30]. Reverse transcriptase real time quantitative PCR (RTqPCR) Changes in gene expression discovered by RNA-seq were previously validated by analysis of a subset of both up and down-regulated genes using RTqPCR [21,25]. Here we analyzed the mRNA levels of a subset of cell cycle control genes using RTqPCR following previously published methods [25]. Oligonucleotide primers used for RTqPCR are given in Table 1. 10.1371/journal.pone.0287858.t001 Table 1 Oligonucleotide primers used for RTqPCR and chromatin immunoprecipitation (ChIP) assay. RTqPCR Gene Primer sequence (5’ → 3’) ccnb1.2 Fwd: GAGGATGCACAAGCAGTCAG Rev: TTTGGGAACTGGGTGTTCCT ccna2 Fwd: TTTGACCTTGCTGCTCCAAC Rev: CGCAGGAAAGGATCAGCATC cdk1 Fwd: GGAACGCCCAACAATGAAGT Rev: CAGGTCCAGCCCATCCTTAT e2f1 Fwd: CGCTGACGTTGTGTATGGTT Rev: TTCCTGTAGGCATTCACGGT a-actinin (reference gene) Fwd: GGACAATTATCCTGCGTTTTGC Rev: CCTTCTTTGGCAGATGTTTCTTC ChIP assay Genomic region Primer sequence (5’ → 3’) thrb TRE Fwd: CCCCTATCCTTGTTCGTCCTC Rev: GCGCTGGGCTGTCCT klf9 synergy module (KSM) Fwd: CCGTCCCTTCTTTTGTGTACATT Rev: GCTGTTCGTGCCACTTTGC thibz TRE Fwd: GGACGCACTAGGGTTAAGTAAGG Rev: TCTCCCAACCCTACAGAGTTCAA ifabp promoter Fwd: CCCTACATTGGTTGAGCCAGTTTT Rev: TCAAAGGCCATGGTGATTGGT thrb exon 5 Fwd: CCCCGAAAGTGAAACTCTAACTCT Rev: CCACACCGAGTCCTCCATTTT Thyroid hormone receptor (TR) chromatin immunoprecipitation (ChIP) and TR ChIP sequencing (TR ChIP-seq) We chose metamorphic climax (NF62) to conduct TR ChIP-seq on tadpole brain because this stage corresponds to the highest circulating plasma TH concentration and thyroid hormone receptor alpha (thra) and thyroid hormone receptor beta (thrb) mRNA levels [1,31]. Also, this stage is when the expression of most regulated genes is highest or lowest, and the rate of tissue transformation is maximal, reflecting active transcriptional regulation by liganded TRs [1]. We isolated chromatin from whole brain (NF62; 5 brains pooled per replicate) and conducted targeted TR ChIP and TR ChIP-seq following our previously published methods [32]. Briefly, we cross-linked the chromatin with formaldehyde, fragmented it by sonication, then precipitated it (5 μg per reaction) using a polyclonal antiserum to Xenopus TRs (PB antiserum provided by Yun-Bo Shi) and the Magna ChIP A/G kit (Millipore) following the manufacturer’s protocol. We analyzed eight replicate ChIP DNA samples using qPCR assays targeting three positive control regions (the Krüppel-like factor 9—klf9 synergy module; the thrb TH response element–TRE; and the thyroid hormone induced bZip protein–thibz—TRE) and two negative control regions (the intestinal fatty acid binding protein [ifabp] promoter and thrb exon 5), then we selected and pooled three replicates with the highest signal/noise ratio. We submitted the sample to the University of Michigan DNA sequencing core along with the input sample for library preparation using the SMARTer ThruPLEX DNA-Seq Kit (Takara), and sequencing in a single lane of an Illumina 4000 Hi-Seq machine at the University of Michigan DNA Sequencing Core. We conducted quality processing on raw reads (fastq) files from each sample to remove poor quality bases and adapter sequences using fastp (v0.20.1). Quality assessments on raw and trimmed reads were performed using FastQC (v0.11.9). We aligned processed reads to the X. tropicalis genome (v9.1), Ensembl release 101, using bwamem (v0.7.17) and default parameters. Alignments were then sorted, indexed using samtools (v1.13), and duplicates marked using Picard (v2.18.2). We filtered unmapped reads, secondary alignments, PCR or optical duplicates, and mappings with a quality below 30 with samtools (v1.13). We created read depth-normalized bigwig files from sorted files using bamCoverage from deepTools (v3.5.1). We identified areas of enriched TR binding using MACS2 (v2.2.6) with the sharp peak type. The regions of the X. tropicalis genome where TR associates in chromatin were visualized using the integrative genome viewer (IGV). For binding sites within 10 kb upstream of a transcription start site (TSS) or within a gene-body, we used ChIPPeakAnno (v3.26.4) to annotate peaks with the corresponding gene name; the distance from the peak to the gene; and the relative position of the peak relative to the gene (ie., upstream vs. internal). We then plotted the genomic distribution of TR binding sites across genomic elements using ChIPPeakAnno (v3.26.4) and ChIPseeker (v1.28.3). We used ChIPseeker (v1.28.3) to create scaled (3 kb upstream and downstream) profile plots of the TR binding sites relative to the TSS and gene-body. Differentially expressed genes (DEGs) were annotated with TR binding site information, defined as enriched peaks with a p-value ≤0.0005, if the peaks were within 10 kb upstream of the gene TSS or within the gene-body. Output from different quality control tools and the peak calling were summarized with MultiQC (v1.11). Select parts of the TR ChIP-seq dataset were previously reported, which included validations using targeted ChIPqPCR assays using chromatin preparations distinct from the chromatin used for ChIP-seq at six known TH-TR target genes (thrb, thibz, klf9, growth arrest and DNA damage inducible gamma—gadd45g, ten eleven translocase 2—tet2, and tet3) and two negative control regions (ifabp promoter, thrb exon 5) [32]. We also conducted targeted ChIPqPCR assays to validate TR ChIP-seq peaks at three uncharacterized loci (S2 Fig). The TR ChIP-seq dataset has been deposited in GEO (GSE214697). Statistical analysis We analyzed RTqPCR data using SYSTAT (v. 13; Systat Software, San Jose, CA). We used one-way ANOVA followed by Fisher’s least significant difference (Fisher’s LSD) post hoc test (a = 0.05). Derived values were Log10-transformed before statistical analysis if the variances were found to be heterogeneous. Results RNA-sequencing analysis in tadpole brain during spontaneous metamorphosis We conducted pairwise analysis of gene expression changes in tadpole brain for five developmental intervals: premetamorphosis to prometamorphosis (NF50-56; i.e., early prometamorphosis); prometamorphosis to metamorphic climax (NF56-62; i.e., late prometamorphosis); metamorphic climax to the completion of metamorphosis (NF62-66); premetamorphosis to metamorphic climax (NF50-62); and premetamorphosis to the completion of metamorphosis (NF50-66) (Tables 2 and S1). There were 5561 unique genes whose mRNA level changed in one or more of the developmental interval comparisons, which represents 26% of the protein coding genes in X. tropicalis (based on an estimated 21,000 genes) [33]. The numbers of differentially expressed genes (DEGs) with pairwise developmental stage comparisons are given in Fig 1A and Tables 2 and S1. We plotted the RNA-seq counts for four known T3-regulated genes and found the expected changes during metamorphosis (Fig 1B). Heatmaps depicting expression changes during metamorphosis for the top 100 upregulated and the top 100 downregulated genes (determined during the developmental interval NF50-NF62) are shown in Fig 1C. 10.1371/journal.pone.0287858.g001 Fig 1 RNA-sequencing analysis in tadpole brain during metamorphosis. A. Pie charts with pairwise comparisons of differentially expressed genes (DEGs) between four stages of metamorphosis. Numbers in red areas of the pie charts represent genes upregulated in Stage B, while the numbers in the blue areas of the pie charts represent genes downregulated in Stage B. Gene regulation changes increased as metamorphosis progressed, with roughly half the genes upregulated and half downregulated. Xenopus illustrations © Natalya Zahn (2022), source Xenbase (www.xenbase.org RRID:SCR_003280) [34]. B. Expression patterns of four known TH-regulated genes in tadpole brain (preoptic area/thalamus/hypothalamus) analyzed at four Nieuwkoop-Faber (NF) stages of spontaneous metamorphosis. Plotted are data from an RNA-seq experiment (n = 3/NF developmental stage). thrb–thyroid hormone receptor b; dio3 –monodiodinase type 3; thibz–thyroid hormone induced bZip protein; klf9 –Krüppel-like factor 9. All except dio3 have been shown to be directly regulated by TH-TR [35–39], which is also supported by the current TR ChIP-seq experiment. C. Heatmaps showing expression changes for the top 100 upregulated and the top 100 downregulated genes (determined during the developmental interval NF50-NF62; S1 Table). 10.1371/journal.pone.0287858.t002 Table 2 Top twenty genes with annotation upregulated from premetamorphosis (NF50) to metamorphic climax (NF62) in tadpole brain. Gene symbol Gene title ENSEMBL Log2fold change* Adjusted p value dio3 deiodinase, iodothyronine, type 3 ENSXETG00000036613 3.275125254 2.94E-40 avt arginine vasotocin** ENSXETG00000017979 2.328280168 1.28E-20 thibz thyroid hormone induced bZip protein ENSXETG00000007922 2.296353968 1.19E-09 spock2 indolethylamine N-methyltransferase ENSXETG00000026030 2.28515105 4.84E-12 vstm5 V-set and transmembrane domain containing 5 ENSXETG00000017117 2.266298133 5.32E-29 phf21a BRAF35-HDAC complex protein BHC80 ENSXETG00000011757 2.214433985 2.69E-32 inhbb inhibin subunit beta B ENSXETG00000036045 2.200045605 1.99E-09 pmepa1 prostate transmembrane protein, androgen induced 1 ENSXETG00000011263 2.195093687 1.46E-24 apoe apolipoprotein E ENSXETG00000035511 2.127879577 8.93E-31 dao D-amino-acid oxidase ENSXETG00000022308 2.112579694 3.37E-24 cyp27a1 cytochrome P450 family 27 subfamily A member 1 ENSXETG00000025024 2.083714291 1.02E-08 p2rx5 purinergic receptor P2X, ligand gated ion channel, 5 ENSXETG00000032944 1.98969067 2.89E-10 cela1.2 chymotrypsin like elastase 1, gene 2 ENSXETG00000039717 1.976354459 2.54E-07 nlrx1 NLR family member X1 ENSXETG00000021586 1.949234863 5.22E-36 arhgap24 Rho GTPase activating protein 24 ENSXETG00000006200 1.887018231 7.6E-09 bhlhe41 basic helix-loop-helix family member e41 ENSXETG00000023482 1.884565228 5.58E-17 arhgap45 Rho GTPase activating protein 45 ENSXETG00000011260 1.877904 6.66E-22 chrm2 cholinergic eceptor muscarinic 2 ENSXETG00000030443 1.877869602 1.87E-12 sbf2 zinc finger protein 816 ENSXETG00000039440 1.825679583 4.21E-11 trpm8 transient receptor potential cation channel, subfamily M, member 8 ENSXETG00000026873 1.824899962 3.94E-15 * Log2fold change of 1 = 2 fold increase. ** This gene is mislabeled in the genome database as arginine vasopressin (avp) which is a mammalian gene. The amphibian gene is arginine vasotocin (avt). The top 20 genes with annotation that were upregulated from NF50 to NF62 are listed in Table 2, and the top 20 genes downregulated during this interval are listed in Table 3. Not included in Table 2 are two hemoglobin subunit genes (hba1 and hbg1) that were strongly upregulated during metamorphosis, and are likely derived from residual blood remaining in the brain at the time of harvest. During amphibian metamorphosis there is a switch from larval to adult erythrocytes that express different globin genes, with hba1 and hbg1 expressed in adult but not in larval cells [40,41]. 10.1371/journal.pone.0287858.t003 Table 3 Top twenty genes with annotation downregulated from premetamorphosis (NF50) to metamorphic climax (NF62) in tadpole brain. Gene symbol Gene title ENSEMBL Log2fold change* Adjusted p value pla2g2e phospholipase A2 group IIE ENSXETG00000037667 -4.444463949 3.8805E-38 col6a6 ribonucleotide reductase M2, gene 2 ENSXETG00000036052 -4.106238088 5.074E-114 crisp1.7 cysteine-rich secretory protein 1 gene 7 ENSXETG00000011152 -3.986829953 1.9676E-29 cdca7 cell division cycle associated 7 ENSXETG00000007518 -3.541269862 1.1933E-55 cct6a component of the chaperonin-containing T-complex (TRiC) ENSXETG00000006636 -3.498033622 1.8951E-53 sapcd1 suppressor APC domain containing 1 ENSXETG00000041061 -3.463026967 4.1714E-66 eef1a1o eukaryotic translation elongation factor 1 alpha 1, oocyte form ENSXETG00000022604 -3.330857635 1.5064E-59 neurog1 neurogenin 1 ENSXETG00000007898 -3.307301672 1.6585E-31 smc2 structural maintenance of chromosomes 2 ENSXETG00000010209 -3.307001586 1.9573E-74 pbk PDZ binding kinase ENSXETG00000020438 -3.236977659 2.3562E-47 kif20a kinesin family member 20A ENSXETG00000006721 -3.193016724 2.1544E-40 cdk2 cyclin-dependent kinase 2 ENSXETG00000009444 -3.185053334 1.8627E-44 cdk1 cyclin-dependent kinase 1 ENSXETG00000003123 -3.183432448 4.228E-56 nusap1 nucleolar and spindle associated protein 1 ENSXETG00000027499 -3.154574577 3.32E-68 kif23 kinesin family member 23 ENSXETG00000007829 -3.15195982 2.6852E-47 cdca8 cell division cycle associated 8 ENSXETG00000002167 -3.1382786 4.1909E-31 ccna2 cyclin A2 ENSXETG00000001016 -3.121308956 6.1124E-44 dlc putative ortholog of delta-like protein C precursor ENSXETG00000002875 -3.090104341 1.629E-79 atad2 ATPase family, AAA domain containing 2 ENSXETG00000023215 -3.068914576 3.5354E-28 hmmr hyaluronan-mediated motility receptor (RHAMM) ENSXETG00000007906 -3.065327172 5.3941E-46 * Log2fold change of -1 = 50% decrease. Patterns of gene regulation during metamorphosis Clustering analysis revealed five major patterns (C1 –C5) of gene regulation during metamorphosis (Fig 2A). Two clusters (C1 and C2) include genes that were downregulated during metamorphosis. The C1 genes, which represent 34.4% of all genes regulated during metamorphosis, were downregulated from NF50 to NF62, but trended up from NF62 to NF66. The C2 genes, which represent 5.5% of all regulated genes, were downregulated from NF50 to NF62 and remained low at NF66. Two clusters (C4 and C5) include genes that were upregulated during metamorphosis. The C4 genes, which represent 34.1% of all genes regulated during metamorphosis, were upregulated from NF50 to NF62, but trended down slightly from NF62 to NF66. The C5 genes, which represent 17.1% of all genes regulated during metamorphosis, were upregulated from NF50 to NF62 and remained elevated at NF66. The C3 genes, which represent 8.9% of all genes regulated during metamorphosis, include genes that were strongly upregulated from NF50 to NF56, remained elevated NF56 to NF62, but then were downregulated from NF62 to NF66. 10.1371/journal.pone.0287858.g002 Fig 2 Patterns of gene regulation in X. tropicalis tadpole brain during spontaneous metamorphosis or after T3 treatment analyzed by RNA-sequencing. A. Clustering analysis showing five patterns of gene expression changes across four stages of metamorphosis. B. Venn diagram showing numbers of genes regulated with overlaps in three developmental stage comparisons. C. Venn diagram showing the overlap between all genes that changed expression during spontaneous metamorphosis with genes induced or repressed by exogenous T3 in premetamorphic tadpoles. Pair-wise comparisons of gene expression changes at four stages of metamorphosis These data are summarized in Table 4. Early prometamorphosis (NF50-56): We identified 1573 genes whose mRNAs levels changed during this interval, with gene downregulation predominating (32% increased, 68% decreased). This developmental period represents the initiation of metamorphosis, and is characterized by rising plasma TH titers, and external changes like the growth and differentiation of the hindlimbs. Late prometamorphosis (NF56-62): We identified 1713 genes whose mRNAs levels changed during this interval, with similar numbers up- and downregulated (53% increased, 47% decreased). This developmental period is characterized by sharply increasing plasma TH titers that peak at metamorphic climax, and external changes like continued hindlimb growth and differentiation, forelimb emergence and growth, intestinal remodeling, cranial re-structuring and the beginning of gill and tail resorption. Metamorphic climax to the completion of metamorphosis (NF62-66): We identified 364 genes whose mRNA levels changed during this interval, 46% up- and 54% downregulated. This developmental period is characterized by a sharp decline in plasma TH titers, completion of resorption of the tail and gills, cranial restructuring, and the emergence of the juvenile frog. Premetamorphosis to metamorphic climax (NF50-62): During the entire period encompassing premetamorphosis to the climax of metamorphosis there were 4005 genes whose mRNA levels changed, 45% up- and 55% downregulated. This developmental period is characterized by a large increase in plasma TH titers from nondetectable to peak concentration, and the progression of the morphological changes that characterize the transformation of the premetamorphic tadpole into the juvenile frog. Premetamorphosis to the completion of metamorphosis (NF50-66): We identified 4185 genes whose mRNA levels changed during the entire interval from the beginning to the end of metamorphosis, with 46% up- and 54% downregulated. 10.1371/journal.pone.0287858.t004 Table 4 Comparisons of gene counts from RNA-sequencing experiments. Numbers of differentially expressed genes Comparison All Up Down Premetamorphosis (NF50) vs prometamorphosis (NF56) 1573 507 1066 Prometamorphosis (NF56) vs metamorphic climax (NF62) 1713 907 806 Metamorphic climax (NF62) vs metamorphic frog (NF66) 364 168 196 Prometamorphosis (NF56) vs metamorphic frog (NF66) 1953 1043 910 Premetamorphosis (NF50) vs metamorphic climax (NF62) 4005 1812 2193 Premetamorphosis (NF50) vs metamorphic frog (NF66) 4185 1923 2262 Number of unique genes regulated during metamorphosis* 5561 T3 treatment of premetamorphic (NF50) tadpoles 2678 1396 1282 Number of unique genes whose mRNA levels change during metamorphosis plus after T3 treatment of premetamorphic tadpoles 6957 Number of genes whose mRNA levels change during metamorphosis and after T3 treatment of premetamorphic tadpoles 1269 Number of genes whose mRNA levels are regulated by T3 in premetamorphic tadpoles but do not change during metamorphosis 1409 Number of genes whose mRNA levels change during metamorphosis but are not regulated by T3 in premetamorphic tadpoles 4297 *These are the unique genes that change their mRNA levels in one or more of the pair-wise developmental interval comparisons. See also S1 and S2 Tables for a full accounting of the RNA-seq data. Overlap: We also looked at the numbers of DEGs that were common to pairs of selected developmental intervals. Comparing the developmental intervals NF50-56 and NF56-62, there were 520 common DEGs, which represented 33.1% and 30.4%, respectively, of all genes regulated during these intervals (Fig 2B and S2 Table). Thus, about one third of the genes whose mRNA levels changed during early prometamorphosis also showed changes during late prometamorphosis/climax; i.e., the mRNA levels for these genes continued to increase or decrease as metamorphosis proceeded to climax. Comparing the developmental intervals NF56-62 and NF62-66, there were 171 common DEGs, which represented 10% and 47%, respectively, of all genes regulated during these intervals (Fig 2B). Thus, only one tenth of the genes whose mRNA levels changed during late prometamorphosis/climax also showed changes from climax to the completion of metamorphosis. Comparing the developmental intervals NF50-56 and NF62-66, there were 112 common DEGs, which represented 7.1% and 30.8%, respectively, of all genes regulated during these intervals (Fig 2B). There were 49 genes that were common between the three developmental intervals. These comparisons show that the expression level of most genes at metamorphic climax, high or low, is maintained in the juvenile frog. RNA-sequencing analysis in brain of premetamorphic tadpoles treated with T3 Treatment of premetamorphic tadpoles (NF54) with T3 (5 nM) for 16 h caused statistically significant changes in the mRNA levels of 2678 genes in tadpole brain (preoptic area/thalamus/hypothalamus); 1396 (52%) were induced and 1282 (48%) were repressed (Fig 2C and Tables 2 and S1). This number represents 13% of the protein coding genes in X. tropicalis. These numbers differ slightly from our previously published estimates [21] owing to our use of the latest genome build in the current study. Overlap of the developmental and T3-induced gene regulation programs Of the unique genes whose mRNA levels changed during the entire metamorphic period (NF50 to NF66; 5561 genes), 1269 genes (22.82%) were common with genes regulated by exogenous T3 in premetamorphic tadpole brain (Fig 2C and Tables 2 and S3). This percentage was similar (22.5%) when comparing the interval NF50-56 when circulating TH concentration is low, but greater (895 of 4005 genes = 33.42%) for the interval NF50-62 when circulating TH concentration increases and peaks (Fig 2B and Tables 2 and S2). Thus, 77% of the genes expressed in tadpole brain whose mRNA levels changed during the entire metamorphic period, and 66% of the genes that changed during NF50-62 were not modulated by exogenous T3 in premetamorphic tadpole brain. On the other hand, of the genes that were regulated by exogenous T3 in premetamorphic tadpole brain (2678 genes), the mRNA levels of 47.39% of these changed during spontaneous metamorphosis. Thus, over half of the genes that are regulated by exogenous T3 in premetamorphic tadpole brain do not change their mRNA levels during metamorphosis. We also analyzed the direction of change in gene expression, up- or downregulated, for the genes that overlapped between the developmental intervals and the T3-induced genes. The concordance was 49.2% for NF50-56, 81% for NF56-62, and 67.2% for NF50-62 (S3 Table). Thus, depending on the developmental interval analyzed, ∼20% to 50% of the genes that change their expression during spontaneous metamorphosis do so in a direction opposite to that caused by treatment of premetamorphic tadpoles with T3. TR association in chromatin in tadpole brain identified by ChIP-sequencing We identified 6302 unique peaks by TR ChIP-seq analysis of brain chromatin from metamorphic climax stage (NF62) tadpoles (hereafter ‘TR peaks’; Tables 5 and S4). The peaks were broadly distributed across each of the 10 X. tropicalis chromosomes (Fig 3A). Our prior validation of the TR ChIP-seq data using targeted ChIP assay confirmed TR association in chromatin at the TRE regions of thrb, thibz, klf9, gadd45g, tet2 and tet3; conversely, there were no TR peaks at two negative control regions (ifabp promoter and thrb exon 5) [32]. Notably, the immediate early gene klf9, whose mRNA level increases at the onset of metamorphosis and shows rapid induction kinetics in response to exogenous T3 in premetamorphic tadpoles [38,42], exhibited a large TR peak at the previously characterized and evolutionarily conserved TRE located in an upstream enhancer (the klf9 synergy module–KSM—located ∼6 kb upstream of the TSS; S3 Fig) [39]. We also discovered additional, previously uncharacterized TR peaks at klf9: one farther upstream to the TSS (∼ -7kb), one just upstream of the TSS, one overlapping the TSS, and one within the gene body (S3 Fig and S5 Table). 10.1371/journal.pone.0287858.g003 Fig 3 Distribution across the genome of thyroid hormone receptor ChIP-seq peaks in X. tropicalis neural cells. A. TR ChIP-seq peaks were found across the entire genome and uniformly distributed on the 10 chromosomes of X. tropicalis. B. A majority of TR ChIP-seq peaks were located +/- 1 kb from the transcription start sites (TSS) of genes. C. Pie chart showing the distribution of TR ChIP-seq peaks across the X. tropicalis genome by genomic feature. D. IGV genome browser tracks showing the locations of TR ChIP-seq peaks (TR peaks; light blue bars) at four loci. Shown are two genes that were upregulated (nucleus accumbens associated 1—nacc1, peak range 0–291; thyroid hormone induced bZip protein—thibz, peak range 0–250) and 2 that were downregulated (aurora kinase A—aurka, peak range 0–326; E2F transcription factor 5—e2f5, peak range 0–233) during metamorphosis. Gene structures are shown in dark blue below the genome tracks. 10.1371/journal.pone.0287858.t005 Table 5 Accounting of TR ChIP-seq peaks. TR ChIP-seq dataset Comparison Number Total number of TR peaks 6302 Number of TR peaks associated with genes 5098 Number of genes* with associated TR peaks 6512 Number of genes without TR peaks 13783 Number of TR peaks with two or more genes 1546 Number of genes with multiple TR peaks 921 Number TR peaks not associated genes 1336 Percentage of TR peaks associated with genes 80.89% TR ChIP-seq peaks at differentially expressed genes (DEGs) and non-DEGs Number of DEGs from beginning to the end of metamorphosis with TR peaks 1355 Number of DEGs from premetamorphosis (NF50) to metamorphic climax (NF62) with multiple TR peaks 207 Number of DEGs from beginning to the end of metamorphosis without TR peaks 4206 Number of non-DEGs from beginning to the end of metamorphosis with TR peaks 3743 Number of DEGs after T3 treatment of premetamorphic (NF50) tadpoles with TR peaks 703 Number of DEGs after T3 treatment of premetamorphic (NF50) tadpoles without TR peaks 1975 * Some TR peaks are associated with more than one gene, and some genes have mulitple TR peaks See S4, S5 and S6 Tables for a full accounting of the TR ChIP-seq data. By setting the window of discovery 10 kb upstream of TSSs and within gene bodies, we found a total of 5098 TR peaks associated with unique genes (which represented 31.29% of the 18881 genes included in the dataset; Tables 5 and S5). Thus, most of the total unique TR peaks (5098 of 6302; 80.89%) were associated with genes using this discovery window. Some TR peaks were associated with more than one gene (1546 peaks), and some genes had multiple TR peaks (921 genes; Tables 5 and S5 and S6). A majority of the TR peaks were clustered +/- 1 kb from the TSSs of the associated genes (Fig 3B); 62.92% were found 0–3 kb upstream, 55.66% 4 k DE genes. This represents > 20% of the total number of genes. In fact, this is so high that ANOVA (and more generally parametric models of variance) fail and should not be used. This increases the level of false positives and negatives. If the comparison was indeed NF50 vs NF66, this is a mistake. But if the list of 4k DE genes is the union of two comparisons (untreated vs T3 treatment NF50, and unterated vs T3 treated NF66), this is a different story. This would raise additional questions, but it is at least technologically sound. This is more than an important technical point, and unfortunately, it is not addressed at all in the manuscript. - The experimental quality of the data set produced sounds reasonable, but the authors fail to provide additional supporting experimental to help the reader estimate its quality. Not that the work hasn't been done correctly, but ChIP-Seq is notoriously difficult to run, and there are numerous (and unavoidable) caveats at the technical and analysis levels. As the manuscript stands, there are very little data to help the reader estimate how far the data the support the conclusions reached by the authors. This is particularly true for the ChIPSeq, which are always very noisy and where peak height is close to background level. - Even if this type of comparisons is popular, the comparison of the ChipSeq peaks with the DE genesis highly speculative. The simple fact that a peak is located in a gene (DE or not) is not evidence of direct regulation. The discussion about this point is long, and they rightfully cite the seminal work of Fullwood et al. So why do the authors label this section "identification of direct thyroid hormone target genes using TR Chrip-seq" ? Can we really say this, without additional experimental data (at least 3C or 4C) ? Minor points - The tool CLUSTRX looks interesting, and the visual output quite telling. Its source should be referenced. Reviewer #2: The study by Raj and colleagues is the first genome-wide analysis of gene expression changes for any tadpole tissue spanning the start to end of natural metamorphosis. Theirs is also the first study to directly compare these changes in spontaneous metamorphosis with those induced precociously by TH, a commonly-used approach. A strong aspect of the current study is that the main conclusions derived from this genome-wide analysis are supported by findings from previous smaller-scale studies. For example, one of the more interesting findings of the study was that, like previous smaller-scale studies, most DEGs in tadpole brain during spontaneous metamorphosis or after T3 treatment did not have TR peaks. The authors put forth several quite plausible explanations for these observations. Another very interesting observation also supported by previous smaller-scale studies was that two thirds of the genes that changed their mRNA levels during spontaneous metamorphosis were not modulated by exogenous T3 in premetamorphic tadpole brain. The authors present several plausible technical and biological explanations for these findings, and conclude that “experiments conducted with exogenous T3 may not be accurate representations of normal developmental processes.” While it is somewhat disappointing that current methods of inducing metamorphosis with T3 and evaluating gene expression do not appear to correspond as well as might be expected with spontaneous metamorphosis, this information is in fact extremely useful to other researchers who use on this method to draw (possibly erroneous) conclusions about natural metamorphosis. One might wonder if inducing metamorphosis instead with a combination of T3 + glucocorticoid may resolve this intriguing discrepancy? Reviewer #3: In this manuscript Raj et al., have done three complementary experiments on Xenopus metamorphosis: (i) they analysed the gene expression changes in the central nervous system at 4 different time point during spontaneous metamorphosis ; (ii) they perform a classical T3 treatment at 5nM during premetamorphosis and look the gene expression change after 16hr; (iii) last they performed a ChIP-seq experiment to detect the regions of chromatin in which the TRs are binding. These three experiments have been on the central nervous system. The main surprise come from the low number of genes that are found in the three methods: As mentioned by the authors, only half of the genes modulated by treatment of premetamorphic tadpoles with TH changed expression during metamorphosis and these represented only 33% of the genes whose mRNA levels changed during metamorphosis. Similarly, only 24% of genes whose mRNA levels changed during metamorphosis had TR ChIP-seq peaks suggesting that they are direct targets. In the Discussion the authors list all the factors that can explain these low numbers. This is a very interesting paper, whose conclusions will probably be disliked by many people working on this field, but it is very useful because it clearly reveals the limits of our assays. Despite our fantastic technical abilities, we are still doing very crude and brutal experiments and by doing that we strongly interfere with the systems we are studying. In addition to this very strong interest, I found the paper well written and the experiments convincing. My main problem is in fact the figures and tables that I find really out of step in terms of quality compared to the text of the paper. Figure 1B and the relevant text lines 217-228 are a bizarre way of presenting the simple fact that the authors have analyzed 4 stages: I would rather present the developmental serie and would link the various stages to show the comparisons that have been made. On Figure 2 and 3 I would add heat maps of the 100 or 200 of top genes and would find a way to visualize the genes that are common between spontaneous metamorphosis and TH regulated. Simply put those genes in color in each of the two heat maps. Figure 4, 5 and 6 are quite boring images of individual genes. I would rather suggest presenting with a scheme the various pathways that are discussed (e.g. cell cycle regulation) and then, again to visualize on each of them the genes that are found by the three technics. There is honestly a strong effort of illustration to do in order to sustain the very nice findings and the great amount of novel information that is present in this paper. Remember that a good part of bioinformatic is to find ways to illustrate very complex and extensive set of information. Also, there is in my point of view one figure missing. At a minimum we need a Vent diagram to compare the number (and %) of genes that are found associated by each of the three methods. These numbers are discussed lines 294-308 and 400-412 and in the Discussion but we need to see in one clear figure how they relate to each other. What is the number of genes that are regulated/associated in the 3 experiments? Are they TH signalling genes? Could we detect that way new genes that would be very important but have escape attention ? Minor points There is an inversion of Figure 5 and 6 with their legends. Lines 241-244 and Figure 2C: there is not a single gene differentially regulated (positively or negatively) in the three stages? This is curious… Fig 2D I would be curious to have few examples of the most interesting genes in each of the C1 tp C5 categories. What I have seen in other paper is the line of 3-4 key gene indicated in color with their names. That would be useful. In our case (in fish) we found many genes implicated in metabolism and a major shift between glycolysis and TCA cycle at different stages during metamorphosis. I think the analysis is a bit short on this point, especially given the importance of the brain as an organ for metabolism also. The link with the appetite control and adipokynes would be super interesting to discuss. ********** 6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files. If you choose “no”, your identity will remain anonymous but your review may still be made public. Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy. Reviewer #1: No Reviewer #2: Yes: Alexander M Schreiber Reviewer #3: No ********** [NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.] While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step. 10.1371/journal.pone.0287858.r002 Author response to Decision Letter 0 Submission Version1 6 Jun 2023 Responses to Reviewer Comments We thank the reviewers for their careful reading of our manuscript and their constructive comments. We have addressed each of the points that were raised as described below. Changes that we made in response to reviewer comments are highlighted in yellow in the revised manuscript file. We also revised the figures as suggested and we rearranged parts of the manuscript to accommodate these changes (we shifted some figure legends and tables). Reviewer #1 The manuscript proposed by Raj, Sifuentes, Kyono and Denver proposes to address an important question of thyroid hormones action in brain during metamorphosis, in the amphibian model Xenopus tropicalis. By combining RNASeq and ChIPseq data, the authors aim to compare transcriptional variations during natural metamorphosis versus following TH treatment. Even if the quality of the writing is certain, the manuscript is far from reaching publication quality. Several subsections of the results section are redundant with other parts of the manuscript. Question: The results description is very concise. Too concise. For example, the comparison of NF50 vs NF66 stages (lines 291-293) reveals > 4 k DE genes. This represents > 20% of the total number of genes. In fact, this is so high that ANOVA (and more generally parametric models of variance) fail and should not be used. This increases the level of false positives and negatives. If the comparison was indeed NF50 vs NF66, this is a mistake. But if the list of 4k DE genes is the union of two comparisons (untreated vs T3 treatment NF50, and unterated vs T3 treated NF66), this is a different story. This would raise additional questions, but it is at least technologically sound. This is more than an important technical point, and unfortunately, it is not addressed at all in the manuscript. Answer: Please note that we used one-way ANOVA only to analyze the targeted RTqPCR data as we indicated in the Materials and Methods. As described in our original publications (Kyono et al., 2020; Wen et al., 2019), we conducted differential expression analysis using the R package DESeq2 (v1.22.0) to identify differentially expressed genes. We reanalyzed the data for the current manuscript using the most recent Xenopus tropicalis genome build (v.9.1). The list of 4K DEGs is not a union of two comparisons as the reviewer indicates; that is, it did not involve the T3 treated vs untreated (NF50) dataset, only the developmental series (NF50, NF56, NF62 and NF66). And there was no T3 treatment of NF66 animals, only of NF50. The NF66 stage is the completion of metamorphosis, and T3 treatment has no developmental/morphological effects in the newly metamorphosed frogs.   Question: The experimental quality of the data set produced sounds reasonable, but the authors fail to provide additional supporting experimental to help the reader estimate its quality. Not that the work hasn't been done correctly, but ChIP-Seq is notoriously difficult to run, and there are numerous (and unavoidable) caveats at the technical and analysis levels. As the manuscript stands, there are very little data to help the reader estimate how far the data the support the conclusions reached by the authors. This is particularly true for the ChIPSeq, which are always very noisy and where peak height is close to background level. Answer: We agree that ChIP-seq (and related techniques) are technically challenging. We have a track record for conducting such studies using antibodies (ChIP-seq), chromatin streptavidin precipitation (ChSPseq) and methyl capture sequencing (MethylCap-seq) (Avila-Mendoza et al., 2020; Knoedler et al., 2017; Kyono et al., 2020). We indicated in the manuscript (original lines 205-206, revised lines 195-201) that we provided targeted ChIP-qPCR validations for the TR ChIP-seq experiment in a previous manuscript published in Endocrinology (Raj et al., 2020) (see both the main Endocrinology manuscript and the supporting information for the manuscript). These included six TH-TR target genes (thrb, thibz, klf9, gadd45�, tet2, tet3) and two negative control regions (ifabp promoter, thrb exon 5). We now include validations for three additional genomic regions in the supporting information of the current manuscript (Supplemental Fig. 2). Lastly, we used a high stringency FDR cutoff (0.01) for our TR ChIP-seq analysis which reduced chances for false positives. Question: Even if this type of comparisons is popular, the comparison of the ChipSeq peaks with the DE genesis highly speculative. The simple fact that a peak is located in a gene (DE or not) is not evidence of direct regulation. The discussion about this point is long, and they rightfully cite the seminal work of Fullwood et al. So why do the authors label this section "identification of direct thyroid hormone target genes using TR Chrip-seq" ? Can we really say this, without additional experimental data (at least 3C or 4C) ? Answer: We agree with the reviewer. We think that we are forthright in describing how and why we selected the criteria that we did to assign TR peaks to genes, and we discussed the caveats for assigning putative TREs (TR peaks) to genes (as the reviewer mentioned). We note that all previously validated TREs in Xenopus have been found within genes, within 5’ UTRs or in the 5’ flanking region, which is why we selected the genomic range (10 kb upstream of the TSS and within gene bodies) that we did to provisionally assign peaks to genes. We agree with the reviewer that more direct evidence is needed to assign functional TREs to genes. We are cognizant of such issues, and earlier we used ChIA-PET data that we validated by 3C assay to show that an upstream TRE in the Xenopus Klf9 gene interacts with the transcription start site (Bagamasbad et al., 2015), and we were the first to use CRISPR-Cas9 genome editing to mutate a putative TRE to demonstrate its transactivation function (in the mouse Dnmt3a gene) (Kyono et al., 2016). We agree that this kind of evidence is required to formally assign a regulatory element to a gene, and that it is premature to conclude that genes with proximate TR peaks are direct thyroid hormone target genes. Therefore, we have changed this label and modified some of the text. Minor points Question: The tool CLUSTRX looks interesting, and the visual output quite telling. Its source should be referenced. Answer: The software is called Clust (not CLUSTRX; this was a mistake, thank you). We have now included the reference in the manuscript. Reviewer #2 The study by Raj and colleagues is the first genome-wide analysis of gene expression changes for any tadpole tissue spanning the start to end of natural metamorphosis. Theirs is also the first study to directly compare these changes in spontaneous metamorphosis with those induced precociously by TH, a commonly-used approach. A strong aspect of the current study is that the main conclusions derived from this genome-wide analysis are supported by findings from previous smaller-scale studies. For example, one of the more interesting findings of the study was that, like previous smaller-scale studies, most DEGs in tadpole brain during spontaneous metamorphosis or after T3 treatment did not have TR peaks. The authors put forth several quite plausible explanations for these observations. Another very interesting observation also supported by previous smaller-scale studies was that two thirds of the genes that changed their mRNA levels during spontaneous metamorphosis were not modulated by exogenous T3 in premetamorphic tadpole brain. The authors present several plausible technical and biological explanations for these findings, and conclude that “experiments conducted with exogenous T3 may not be accurate representations of normal developmental processes.” While it is somewhat disappointing that current methods of inducing metamorphosis with T3 and evaluating gene expression do not appear to correspond as well as might be expected with spontaneous metamorphosis, this information is in fact extremely useful to other researchers who use on this method to draw (possibly erroneous) conclusions about natural metamorphosis. One might wonder if inducing metamorphosis instead with a combination of T3 + glucocorticoid may resolve this intriguing discrepancy? Answer: We thank the reviewer for their supportive comments. It might be interesting to analyze the combined T3 + glucocorticoid treatment. However, we believe that the possible artifactual responses are more likely related to the pharmacological effect of exposing premetamorphic tadpoles to a metamorphic climax stage dose of T3 (or even greater, as we found previously that tadpoles concentrate hormone from their environment, and we discuss in the text). A dose response experiment with T3 (and a time course) might resolve this issue. Reviewer #3 In this manuscript Raj et al., have done three complementary experiments on Xenopus metamorphosis: (i) they analysed the gene expression changes in the central nervous system at 4 different time point during spontaneous metamorphosis ; (ii) they perform a classical T3 treatment at 5nM during premetamorphosis and look the gene expression change after 16hr; (iii) last they performed a ChIP-seq experiment to detect the regions of chromatin in which the TRs are binding. These three experiments have been on the central nervous system. The main surprise come from the low number of genes that are found in the three methods: As mentioned by the authors, only half of the genes modulated by treatment of premetamorphic tadpoles with TH changed expression during metamorphosis and these represented only 33% of the genes whose mRNA levels changed during metamorphosis. Similarly, only 24% of genes whose mRNA levels changed during metamorphosis had TR ChIP-seq peaks suggesting that they are direct targets. In the Discussion the authors list all the factors that can explain these low numbers. This is a very interesting paper, whose conclusions will probably be disliked by many people working on this field, but it is very useful because it clearly reveals the limits of our assays. Despite our fantastic technical abilities, we are still doing very crude and brutal experiments and by doing that we strongly interfere with the systems we are studying. In addition to this very strong interest, I found the paper well written and the experiments convincing. My main problem is in fact the figures and tables that I find really out of step in terms of quality compared to the text of the paper. Answer: We thank the reviewer for their supportive comments. We have now revised the figures and tables as described below. Question: Figure 1B and the relevant text lines 217-228 are a bizarre way of presenting the simple fact that the authors have analyzed 4 stages: I would rather present the developmental series and would link the various stages to show the comparisons that have been made. Answer: The purpose of the original figure 1B was to show the comparisons that we made among the RNA-seq data from the 4 developmental stages, but we agree that this attempt was not the best. To aid readers who may not be familiar with the gross morphological changes that occur during metamorphosis of X. tropicalis we now include tadpole images on the new figure 1B. The text on lines 217-228 (now lines 213-218) simply states the comparisons that we made between the mRNA levels determined by RNA-seq at the 4 developmental stages. Question: On Figure 2 and 3 I would add heat maps of the 100 or 200 of top genes and would find a way to visualize the genes that are common between spontaneous metamorphosis and TH regulated. Simply put those genes in color in each of the two heat maps. Answer: We now provide heatmaps that show the top 100 upregulated and the top 100 downregulated genes during metamorphosis in figure 2. Although these heatmaps do not provide additional information beyond what we presented in the clustering analysis, which includes many more genes (Fig. 1), they may help the reader visualize the gene expression changes that occur during metamorphosis of tadpole brain. However, we cannot generate a heat map comparing expression levels between the developmental series and the +/- T3 treatment since these were two separate experiments. Question: Figure 4, 5 and 6 are quite boring images of individual genes. I would rather suggest presenting with a scheme the various pathways that are discussed (e.g. cell cycle regulation) and then, again to visualize on each of them the genes that are found by the three technics. There is honestly a strong effort of illustration to do in order to sustain the very nice findings and the great amount of novel information that is present in this paper. Remember that a good part of bioinformatic is to find ways to illustrate very complex and extensive set of information. Answer: We now include a KEGG pathway analysis diagram for cell cycle in figure 4, and for ribosome biogenesis in figure 5. We chose to not provide similar pathway figures for the T3 experiment because we do not trust that such analyses are reliable representations of normal developmental processes, as we discuss in the manuscript. We do not provide GO or pathway analyses for the TR ChIP-seq experiment because we cannot determine if genes with proximate TR peaks are actually regulated by the T3-TR complex. Indeed, our analysis showed that 2/3 of the TR peaks associated with genes were found proximate to genes that did not change expression during metamorphosis (∼5/6 of TR peaks were not associated with genes regulated by exogenous T3). Therefore, the only dataset that we think is amenable to such analysis is the developmental series, which includes differentially expressed genes that are likely to be involved in the biological processes under investigation. We decided to keep figures 6 and 7 since they depict changes in the expression of key genes involved with neural differentiation and neuroendocrine function, which we think that researchers working on Xenopus will find interesting. Question: Also, there is in my point of view one figure missing. At a minimum we need a Vent diagram to compare the number (and %) of genes that are found associated by each of the three methods. These numbers are discussed lines 294-308 and 400-412 and in the Discussion but we need to see in one clear figure how they relate to each other. What is the number of genes that are regulated/associated in the 3 experiments? Are they TH signalling genes? Could we detect that way new genes that would be very important but have escape attention? Answer: We now provide additional an Venn diagram in figure 2 comparing differentially expressed genes during spontaneous metamorphosis and after T3 treatment of premetamorphic tadpoles. However, we cannot generate Venn diagrams that also include the TR ChIP-seq data because the measurements are different (i.e., DEGs vs. TR peaks). Table 5 summarizes the TR ChIP-seq data. Most of the TR peaks were in proximity to genes that did not change during metamorphosis or after T3 treatment. These genes may not be regulated by T3-TR at all, or they may not be regulated by T3-TR during the developmental stages that we analyzed; or, the TREs at which the TRs associate regulate genes other than those with which they are in close proximity (see Discussion and our answer to reviewer #1 above). There are several genes involved with T3 action that we highlight in the manuscript. It is possible that other genes in our lists play a role in T3 signaling, but we have no way of determining this from our analyses. We provide extensive supplemental tables that we hope will be a valuable resource for future investigations. Minor points Question: There is an inversion of Figure 5 and 6 with their legends. Answer: Thank you for noticing this error which has now been corrected. Question: Lines 241-244 and Figure 2C: there is not a single gene differentially regulated (positively or negatively) in the three stages? This is curious… Answer: Thank you for noticing this omission which has now been corrected. There were 49 genes that were common between the three developmental intervals (NF50-56, NF56-62 and NF62-66). This number is low in part because the number of genes regulated between NF62-66 was low (364). Question: Fig 2D I would be curious to have few examples of the most interesting genes in each of the C1 tp C5 categories. What I have seen in other paper is the line of 3-4 key gene indicated in color with their names. That would be useful. In our case (in fish) we found many genes implicated in metabolism and a major shift between glycolysis and TCA cycle at different stages during metamorphosis. I think the analysis is a bit short on this point, especially given the importance of the brain as an organ for metabolism also. The link with the appetite control and adipokynes would be super interesting to discuss. Answer: The expression patterns of some of these genes are shown in figures 5-7. We chose these genes to plot because of their previously demonstrated (or hypothesized based on studies in other species) importance for Xenopus development. We think that Xenopus researchers will be interested in the developmental profiles of these genes, which we categorized based on their functions in discrete developmental/cellular signaling pathways. We agree that genes involved with energy balance and the regulation of food intake are interesting, which we present briefly. However, one could make a similar case for reproduction, pituitary regulation, among others. We think that these are interesting discussion points, but to do them justice would require that the manuscript be considerably longer than it already is, and so we think that they are beyond the scope of the current manuscript. References Cited Avila-Mendoza, J., Subramani, A., Sifuentes, C.J., Denver, R.J., 2020. Molecular Mechanisms for Kruppel-Like Factor 13 Actions in Hippocampal Neurons. Molecular Neurobiology 57, 3785-3802. Bagamasbad, P., Bonett, R., Sachs, L., Buisine, N., Raj, S., Knoedler, J., Kyono, Y., Ruan, Y., Ruan, X., Denver, R., 2015. Deciphering the regulatory logic of an ancient, ultraconserved nuclear receptor enhancer module. Molecular Endocrinology 29, 856-872. Knoedler, J.R., Subramani, A., Denver, R.J., 2017. The Kruppel-like factor 9 cistrome in mouse hippocampal neurons reveals predominant transcriptional repression via proximal promoter binding. BMC Genomics 18, 299. Kyono, Y., Raj, S., Sifuentes, C.J., Buisine, N., Sachs, L., Denver, R.J., 2020. DNA methylation dynamics underlie metamorphic gene regulation programs in Xenopus tadpole brain. Dev. Biol. 462, 180-196. Kyono, Y., Subramani, A., Ramadoss, P., Hollenberg, A.N., Bonett, R.M., Denver, R.J., 2016. Liganded thyroid hormone receptors transactivate the DNA methyltransferase 3a gene in mouse neuronal cells. Endocrinology 157, 3647-3657. Raj, S., Kyono, Y., Sifuentes, C.J., Arellanes-Licea, E., Subramani, A., Denver, R.J., 2020. Thyroid hormone induces DNA demethylation in Xenopus tadpole brain. Dryad Digital Repository Deposited 1 July, 2020. Wen, L., He, C., Sifuentes, C.J., Denver, R.J., 2019. Thyroid hormone receptor alpha Is required for thyroid hormone-dependent neural cell proliferation during tadpole metamorphosis. Frontiers in Endocrinology 10. Attachment Submitted filename: Response to Reviewers.docx Click here for additional data file. 10.1371/journal.pone.0287858.r003 Decision Letter 1 Coen Laurent Academic Editor © 2023 Laurent Coen 2023 Laurent Coen https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Submission Version1 14 Jun 2023 Metamorphic gene regulation programs in Xenopus tropicalis tadpole brain PONE-D-23-06337R1 Dear Dr. Denver, We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements. Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication. An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org. If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org. Kind regards, Laurent Coen, Ph.D. Academic Editor PLOS ONE Additional Editor Comments (optional): After a careful reading of the corrections made by the authors, I consider that all the points raised by the reviewers have been taken into account. All the modifications and comments asked by the reviewers have been addressed in the revised version of the manuscript, and important changes in the text and in figs presentation have been done. Responses to comments, rewordings of the text and changes to the figures have clarified the paper's weaknesses pointed out in its initial version, and it can therefore be accepted for publication in its revised form, without the need for further consideration by the reviewers. Reviewers' comments: 10.1371/journal.pone.0287858.r004 Acceptance letter Coen Laurent Academic Editor © 2023 Laurent Coen 2023 Laurent Coen https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 19 Jun 2023 PONE-D-23-06337R1 Metamorphic gene regulation programs in Xenopus tropicalis tadpole brain Dear Dr. Denver: I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department. If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org. If we can help with anything else, please email us at plosone@plos.org. Thank you for submitting your work to PLOS ONE and supporting open access. Kind regards, PLOS ONE Editorial Office Staff on behalf of Dr. Laurent Coen Academic Editor PLOS ONE ==== Refs References 1 Buchholz DR . Xenopus metamorphosis as a model to study thyroid hormone receptor function during vertebrate developmental transitions. Mol Cell Endocrinol. 2017;459 (C ):64–70. doi: 10.1016/j.mce.2017.03.020 PubMed PMID: WOS:000423244800009. 28363743 2 Cheng S-Y , Leonard JL , Davis PJ . Molecular aspects of thyroid hormone actions. Endocr Rev. 2010;31 (2 ):139–70. doi: 10.1210/er.2009-0007 .20051527 3 Brown DD , Cai L . Amphibian metamorphosis. Dev Biol. 2007;306 (1 ):20. doi: 10.1016/j.ydbio.2007.03.021 .17449026 4 Denver RJ . Endocrinology of Complex Life Cycles: Amphibians. Balthazart J , Pfaff DW , Joels M , editors 2017. 5 Kikuyama S , Hasunuma I , Okada R . Development of the hypothalamo-hypophyseal system in amphibians with special reference to metamorphosis. Mol Cell Endocrinol. 2021;524 . doi: 10.1016/j.mce.2020.111143 PubMed PMID: WOS:000617765000001. 33385474 6 Wang Z , Brown DD . A gene expression screen. Proc Natl Acad Sci U S A. 1991;88 (24 ):11505–9. doi: 10.1073/pnas.88.24.11505 1722336 7 Buckbinder L , Brown DD . Thyroid hormone-induced gene expression changes in the developing frog limb. J Biol Chem. 1992;262 :11221–7. PubMed PMID: .1464592 8 Kanamori A , Brown DD . Cultured cells as a model for amphibian metamorphosis. Proc Natl Acad Sci U S A. 1993;90 (13 ):6013–7. doi: 10.1073/pnas.90.13.6013 8327476 9 Shi YB , Brown DD . The earliest changes in gene expression in tadpole intestine induced by thyroid hormone. J Biol Chem. 1993;268 (27 ):20312–7. PubMed PMID: .7690754 10 Wang Z , Brown DD . Thyroid hormone-induced gene expression program for amphibian tail resorption. J Biol Chem. 1993;268 (22 ):16270–8. PubMed PMID: .8344914 11 Denver RJ , Pavgi S , Shi YB . Thyroid hormone-dependent gene expression program for Xenopus neural development. J Biol Chem. 1997;272 (13 ):8179–88. PubMed PMID: ISI:A1997WQ63500018. doi: 10.1074/jbc.272.13.8179 9079635 12 Fu LZ , Das B , Matsuura K , Fujimoto K , Heimeier RA , Shi YB . Genome-wide identification of thyroid hormone receptor targets in the remodeling intestine during Xenopus tropicalis metamorphosis. Scientific Reports. 2017;7 . doi: 10.1038/s41598-017-06679-x PubMed PMID: WOS:000406280000007. 28743885 13 Heimeier RA , Das B , Buchholz DR , Fiorentino M , Shi YB . Studies on Xenopus laevis intestine reveal biological pathways underlying vertebrate gut adaptation from embryo to adult. Genome Biol. 2010;11 (5 ):R55. Epub 2010/05/21. doi: 10.1186/gb-2010-11-5-r55 ; PubMed Central PMCID: PMC2898076.20482879 14 Buchholz DR , Heimeier RA , Das B , Washington T , Shi YB . Pairing morphology with gene expression in thyroid hormone-induced intestinal remodeling and identification of a core set of TH-induced genes across tadpole tissues. Dev Biol. 2007;303 (2 ):576–90. PubMed PMID: ISI:000244991500015. doi: 10.1016/j.ydbio.2006.11.037 17214978 15 Das B , Heimeier RA , Buchholz DR , Shi YB . Identification of direct thyroid hormone response genes reveals the earliest gene regulation programs during frog metamorphosis. J Biol Chem. 2009;284 (49 ):34167–78. doi: 10.1074/jbc.M109.066084 PubMed PMID: WOS:000272165200049. 19801647 16 Sun G , Heimeier RA , Fu L , Hasebe T , Das B , Ishizuya-Oka A , et al . Expression profiling of intestinal tissues implicates tissue-specific genes and pathways essential for thyroid hormone-induced adult stem cell development. Endocrinology. 2013;154 (11 ):4396–407. doi: 10.1210/en.2013-1432 PubMed PMID: WOS:000325935200047. 23970787 17 Kulkarni SS , Buchholz D . Beyond synergy: Corticosterone and thyroid hormone have numerous interaction effects on gene regulation in Xenopus tropicalis tadpoles. Endocrinology. 2012;en.2012–1432; doi: 10.1210/en.2012-1432 22315456 18 Kulkarni SS , Buchholz DR . Developmental programs and endocrine disruption in frog metamorphosis: The perspective from microarray analysis. In: Shi YB , editor. Animal Metamorphosis. Current Topics in Developmental Biology. 1032013. p. 329–64. doi: 10.1016/B978-0-12-385979-2.00012-5 23347525 19 Heimeier RA , Das B , Buchholz DR , Shi YB . The Xenoestrogen Bisphenol A Inhibits Postembryonic Vertebrate Development by Antagonizing Gene Regulation by Thyroid Hormone. Endocrinology. 2009;150 (6 ):2964–73. doi: 10.1210/en.2008-1503 PubMed PMID: ISI:000266256700062. 19228888 20 Searcy BT , Beckstrom-Sternberg SM , Beckstrom-Sternberg JS , Stafford P , Schwendiman AL , Soto-Pena J , et al . Thyroid hormone-dependent development in Xenopus laevis: A sensitive screen of thyroid hormone signaling disruption by municipal wastewater treatment plant effluent. General and Comparative Endocrinology. 2012;176 (3 ):481–92. doi: 10.1016/j.ygcen.2011.12.036 PubMed PMID: WOS:000303700900031. 22248444 21 Wen L , He C , Sifuentes CJ , Denver RJ . Thyroid hormone receptor alpha Is required for thyroid hormone-dependent neural cell proliferation during tadpole metamorphosis. Frontiers in Endocrinology. 2019;10 . doi: 10.3389/fendo.2019.00396 PubMed PMID: WOS:000473273200001. 31316462 22 Tanizaki Y , Shibata Y , Zhang HG , Shi YB . Analysis of thyroid hormone receptor alpha-knockout tadpoles reveals that the activation of cell cycle program is involved in thyroid hormone-induced larval epithelial cell death and adult intestinal stem cell development during Xenopus tropicalis metamorphosis. Thyroid. 2021;31 (1 ):128–42. doi: 10.1089/thy.2020.0022 PubMed PMID: WOS:000547646900001. 32515287 23 Cordero-Veliz C , Larrain J , Faunes F . Transcriptome analysis of the response to thyroid hormone in Xenopus neural stem and progenitor cells. Developmental Dynamics.11 . doi: 10.1002/dvdy.535 PubMed PMID: WOS:000855000000001. 36065982 24 Ta AC , Huang LC , McKeown CR , Bestman JE , Van Keuren-Jensen K , Cline HT . Temporal and spatial transcriptomic dynamics across brain development in Xenopus laevis tadpoles. G3-Genes Genomes Genetics. 2022;12 (1 ):14. doi: 10.1093/g3journal/jkab387 PubMed PMID: WOS:000743786100007. 34751375 25 Kyono Y , Raj S , Sifuentes CJ , Buisine N , Sachs L , Denver RJ . DNA methylation dynamics underlie metamorphic gene regulation programs in Xenopus tadpole brain. Dev Biol. 2020;462 (2 ):180–96. doi: 10.1016/j.ydbio.2020.03.013 PubMed PMID: WOS:000535972600006. 32240642 26 Nieuwkoop PD , Faber J . Normal Table of Xenopus laevis (Daudin). New York: Garland Publishing Inc; 1994. 27 Babicki S , Arndt D , Marcu A , Liang YJ , Grant JR , Maciejewski A , et al . Heatmapper: web-enabled heat mapping for all. Nucleic Acids Research. 2016;44 (W1 ):W147–W53. doi: 10.1093/nar/gkw419 PubMed PMID: WOS:000379786800025. 27190236 28 Abu-Jamous B , Kelly S . Clust: automatic extraction of optimal co-expressed gene clusters from gene expression data. Genome Biology. 2018;19 :172. doi: 10.1186/s13059-018-1536-8 30359297 29 Reimand J , Arak T , Vilo J . g:Profiler-a web server for functional interpretation of gene lists (2011 update). Nucleic Acids Research. 2011;39 :W307–W15. doi: 10.1093/nar/gkr378 PubMed PMID: WOS:000292325300050. 21646343 30 Luo WJ , Brouwer C . PATHVIEW WEB: User friendly pathway visualization and data integration. Eur Neuropsychopharmacol. 2019;29 :S963–S. doi: 10.1016/j.euroneuro.2017.08.323 PubMed PMID: WOS:000462156400452. 31 Shi Y-B . Amphibian metamorphosis: from morphology to molecular biology. New York: Wiley-Liss; 2000. xiv, 288 p. p. 32 Raj S , Kyono Y , Sifuentes CJ , Arellanes-Licea ED , Subramani A , Denver RJ . Thyroid Hormone Induces DNA Demethylation in Xenopus Tadpole Brain. Endocrinology. 2020;161 (11 ):18. doi: 10.1210/endocr/bqaa155 PubMed PMID: WOS:000754315400001. 32865566 33 Hellsten U , Harland RM , Gilchrist MJ , Hendrix D , Jurka J , Kapitonov V , et al . The genome of the Western clawed frog Xenopus tropicalis. Science. 2010;328 (5978 ):633–6. doi: 10.1126/science.1183670 PubMed PMID: WOS:000277159800047. 20431018 34 Zahn N , James-Zorn C , Ponferrada VG , Adams DS , Grzymkowski J , Buchholz DR , et al . Normal Table of Xenopus development: a new graphical resource. Development. 2022;149 (14 ):16. doi: 10.1242/dev.200356 PubMed PMID: WOS:000835484900008. 35833709 35 Ranjan M , Wong J , Shi Y-B . Transcriptional repression of Xenopus TRb gene is mediated by a thyroid hormone response element located near the start site. J Biol Chem. 1994;269 (40 ):24699–705.7929143 36 Machuca I , Esslemont G , Fairclough L , Tata JR . Analysis of structure and expression of the Xenopus thyroid hormone receptor-beta gene to explain its autoinduction. Mol Endocrinol. 1995;9 (1 ):96–107. Epub 1995/01/01. PubMed doi: 10.1210/mend.9.1.7760854 .7760854 37 Furlow JD , Brown DD . In vitro and in vivo analysis of the regulation of a transcription factor gene by thyroid hormone during Xenopus laevis metamorphosis. Molecular Endocrinology. 1999;13 (12 ):2076–89. PubMed PMID: ISI:000084093600009. doi: 10.1210/mend.13.12.0383 10598583 38 Furlow JD , Kanamori A . The transcription factor basic transcription element-binding protein 1 is a direct thyroid hormone response gene in the frog Xenopus laevis. Endocrinology. 2002;143 (9 ):3295–305. doi: 10.1210/en.2002-220126 12193541 39 Bagamasbad P , Bonett R , Sachs L , Buisine N , Raj S , Knoedler J , et al . Deciphering the regulatory logic of an ancient, ultraconserved nuclear receptor enhancer module. Molecular Endocrinology. 2015;29 (6 ):856–72. doi: 10.1210/me.2014-1349 25866873 40 Mukhi S , Cai LQ , Brown DD . Gene switching at Xenopus laevis metamorphosis. Dev Biol. 2010;338 (2 ):117–26. doi: 10.1016/j.ydbio.2009.10.041 PubMed PMID: WOS:000274436000001. 19896938 41 Liao Y , Ma LF , Guo QL , Weigao E , Fang X , Yang L , et al . Cell landscape of larval and adult Xenopus laevis at single-cell resolution. Nature Communications. 2022;13 (1 ):15. doi: 10.1038/s41467-022-31949-2 PubMed PMID: WOS:000848744000001. 35879314 42 Hoopfer ED , Huang LY , Denver RJ . Basic transcription element binding protein is a thyroid hormone-regulated transcription factor expressed during metamorphosis in Xenopus laevis. Dev Growth Diff. 2002;44 (5 ):365–81. PubMed PMID: ISI:000178809800002. 43 Tanizaki Y , Zhang H , Shibata Y , Shi YB . Thyroid hormone receptor alpha controls larval intestinal epithelial cell death by regulating the CDK1 pathway. Commun Biol. 2022;5 (1 ):10. doi: 10.1038/s42003-022-03061-0 PubMed PMID: WOS:000752399600001. 35013537 44 Hones GS , Kerp H , Hoppe C , Kowalczyk M , Zwanziger D , Baba HA , et al . Canonical thyroid hormone receptor beta action stimulates hepatocyte proliferation in male mice. Endocrinology. 2022;163 (3 ). doi: 10.1210/endocr/bqac003 PubMed PMID: WOS:000753119900005. 35038735 45 Kress E , Samarut J , Plateroti M . Thyroid hormones and the control of cell proliferation or cell differentiation: Paradox or duality? Mol Cell Endocrinol. 2009;313 (1–2 ):36–49. doi: 10.1016/j.mce.2009.08.028 PubMed PMID: WOS:000271488700005. 19737599 46 Denver RJ , Hu F , Scanlan TS , Furlow JD . Thyroid hormone receptor subtype specificity for hormone-dependent neurogenesis in Xenopus laevis. Dev Biol. 2009;326 (1 ):155–68. doi: 10.1016/j.ydbio.2008.11.005 .19056375 47 Bonett RM , Hoopfer ED , Denver RJ . Molecular mechanisms of corticosteroid synergy with thyroid hormone during tadpole metamorphosis. General and Comparative Endocrinology. 2010;168 (2 ):209–19. doi: 10.1016/j.ygcen.2010.03.014 PubMed PMID: ISI:000281095000007. 20338173 48 Denver RJ . The molecular basis of thyroid hormone-dependent central nervous system remodeling during amphibian metamorphosis1. Comp Biochem Physiol C Pharmacol Toxicol Endocrinol. 1998;119 (3 ):219–28. doi: 10.1016/S0742-8413(98)00011-5 9826995 49 Kikuyama S , Kawamura K , Tanaka S , Yamamoto K . Aspects of amphibian metamorphosis: hormonal control. Int Rev Cytol. 1993;145 :105–48. doi: 10.1016/s0074-7696(08)60426-x 8500980 50 Paul B , Sterner ZR , Buchholz DR , Shi YB , Sachs LM . Thyroid and corticosteroid signaling in amphibian metamorphosis. Cells. 2022;11 (10 ):23. doi: 10.3390/cells11101595 PubMed PMID: WOS:000804272400001. 35626631 51 Chatonnet F , Guyot R , Benoit G , Flamant F . Genome-wide analysis of thyroid hormone receptors shared and specific functions in neural cells. Proc Natl Acad Sci U S A. 2013;110 (8 ):E766–75. doi: 10.1073/pnas.1210626110 ; PubMed Central PMCID: PMC3581916.23382204 52 Ramadoss P , Abraham BJ , Tsai L , Zhou Y , Costa-e-Sousa RH , Ye F , et al . Novel mechanism of positive versus negative regulation by thyroid hormone receptor beta1 (TRbeta1) identified by genome-wide profiling of binding sites in mouse liver. J Biol Chem. 2014;289 (3 ):1313–28. doi: 10.1074/jbc.M113.521450 ; PubMed Central PMCID: PMC3894317.24288132 53 Grøntved L , Waterfall JJ , Kim DW , Baek S , Sung M-H , Zhao L , et al . Transcriptional activation by the thyroid hormone receptor through ligand-dependent receptor recruitment and chromatin remodelling. Nat Commun. 2015;6 :7048. doi: 10.1038/ncomms8048 25916672 54 Ayers S , Switnicki MP , Angajala A , Lammel J , Arumanayagam AS , Webb P . Genome-wide binding patterns of thyroid hormone receptor beta. Plos One. 2014;9 (2 ):14. doi: 10.1371/journal.pone.0081186 PubMed PMID: WOS:000331706700002. 24558356 55 Dong HY , Yauk CL , Rowan-Carroll A , You SH , Zoeller RT , Lambert L , et al . Identification of thyroid hormone receptor binding sites and target genes using ChIP-on-Chip in developing mouse cerebellum. Plos One. 2009;4 (2 ). doi: 10.1371/journal.pone.0004610 PubMed PMID: WOS:000278086700020. 19240802 56 Zekri Y , Guyot R , Flamant F . An Atlas of Thyroid Hormone Receptors’ Target Genes in Mouse Tissues. International Journal of Molecular Sciences. 2022;23 (19 ):13. doi: 10.3390/ijms231911444 PubMed PMID: WOS:000867726300001. 36232747 57 Galouzis CC , Furlong EEM . Regulating specificity in enhancer-promoter communication. Current Opinion in Cell Biology. 2022;75 :11. doi: 10.1016/j.ceb.2022.01.010 PubMed PMID: WOS:000798971300013. 35240372 58 Ozdemir I , Gambetta MC . The role of insulation in patterning gene expression. Genes. 2019;10 (10 ):22. doi: 10.3390/genes10100767 PubMed PMID: WOS:000498397100038. 31569427 59 Herold M , Bartkuhn M , Renkawitz R . CTCF: insights into insulator function during development. Development. 2012;139 (6 ):1045–57. doi: 10.1242/dev.065268 PubMed PMID: WOS:000300640700001. 22354838 60 Li G , Fullwood MJ , Xu H , Mulawadi FH , Velkov S , Vega V , et al . ChIA-PET tool for comprehensive chromatin interaction analysis with paired-end tag sequencing. Genome Biology. 2010;11 (2 ). doi: 10.1186/gb-2010-11-2-r22 PubMed PMID: ISI:000276434300018. 20181287 61 Fullwood MJ , Ruan YJ . ChIP-Based Methods for the Identification of Long-Range Chromatin Interactions. J Cell Biochem. 2009;107 (1 ):30–9. doi: 10.1002/jcb.22116 PubMed PMID: ISI:000265710100005. 19247990 62 Chatonnet F , Flamant F , Morte B . A temporary compendium of thyroid hormone target genes in brain. Biochim Biophys Acta. 2015;1849 (2 ):122–9. doi: 10.1016/j.bbagrm.2014.05.023 .24882357 63 Thompson CK , Cline HT . Thyroid hormone acts locally to increase neurogenesis, neuronal differentiation, and dendritic arbor elaboration in the tadpole visual system. Journal of Neuroscience. 2016;36 (40 ):10356–75. doi: 10.1523/JNEUROSCI.4147-15.2016 PubMed PMID: WOS:000388365000013. 27707971 64 Thuret R , Auger H , Papalopulu N . Analysis of neural progenitors from embryogenesis to juvenile adult in Xenopus laevis reveals biphasic neurogenesis and continuous lengthening of the cell cycle. Biology Open. 2015;4 (12 ):1772–81. doi: 10.1242/bio.013391 PubMed PMID: WOS:000366672900020. 26621828 65 Shtutman M , Zhurinsky J , Simcha I , Albanese C , D’Amico M , Pestell R , et al . The cyclin D1 gene is a target of the beta-catenin/LEF-1 pathway. Proc Natl Acad Sci U S A. 1999;96 (10 ):5522–7. doi: 10.1073/pnas.96.10.5522 PubMed PMID: WOS:000080246500038. 10318916 66 Davidson G , Niehrs C . Emerging links between CDK cell cycle regulators and Wnt signaling. Trends in Cell Biology. 2010;20 (8 ):453–60. doi: 10.1016/j.tcb.2010.05.002 PubMed PMID: WOS:000281098000003. 20627573 67 Miller LD , Park KS , Guo QBM , Alkharouf NW , Malek RL , Lee NH , et al . Silencing of Wnt signaling and activation of multiple metabolic pathways in response to thyroid hormone-stimulated cell proliferation. Mol Cell Biol. 2001;21 (19 ):6626–39. doi: 10.1128/MCB.21.19.6626-6639.2001 PubMed PMID: WOS:000170842100025. 11533250 68 Guigon CJ , Zhao L , Lu CX , Willingham MC , Cheng SY . Regulation of beta-catenin by a novel nongenomic action of thyroid hormone beta receptor. Mol Cell Biol. 2008;28 (14 ):4598–608. doi: 10.1128/MCB.02192-07 PubMed PMID: WOS:000257613400013. 18474620 69 Perez-Juste G , Aranda A . The cyclin-dependent kinase inhibitor p27(Kip1) is involved in thyroid hormone-mediated neuronal differentiation. J Biol Chem. 1999;274 (8 ):5026–31. PubMed PMID: ISI:000078698200073. doi: 10.1074/jbc.274.8.5026 9988748 70 Vlaicu SI , Tatomir A , Anselmo F , Boodhoo D , Chira R , Rus V , et al . RGC-32 and diseases: the first 20 years. Immunol Res. 2019;67 (2–3 ):267–79. doi: 10.1007/s12026-019-09080-0 PubMed PMID: WOS:000484525000012. 31250246 71 Bender MC , Hu C , Pelletier C , Denver RJ . To eat or not to eat: ontogeny of hypothalamic feeding controls and a role for leptin in modulating life-history transition in amphibian tadpoles. Proceedings of the Royal Society B-Biological Sciences. 2018;285 (1875 ). doi: 10.1098/rspb.2017.2784 PubMed PMID: WOS:000428940600010. 29593109 72 Chen M , Puschmann TB , Marasek P , Inagaki M , Pekna M , Wilhelmsson U , et al . Increased neuronal differentiation of neural progenitor cells derived from phosphovimentin-deficient mice. Molecular Neurobiology. 2018;55 (7 ):5478–89. doi: 10.1007/s12035-017-0759-0 PubMed PMID: WOS:000434805100006. 28956310 73 Wilhelmsson U , Lebkuechner I , Leke R , Marasek P , Yang XG , Antfolk D , et al . Nestin regulates neurogenesis in mice through notch signaling from astrocytes to neural stem cells. Cereb Cortex. 2019;29 (10 ):4050–66. doi: 10.1093/cercor/bhy284 PubMed PMID: WOS:000493326300003. 30605503 74 Kobayashi T, Kageyama R. Expression dynamics and functions of hes factors in development and diseases. In: Taneja R, editor. Current Topics in Developmental Biology. BLH transcription factors in development and disease. 110. San Diego: Elsevier Academic Press Inc; 2014. p. 263–83. 75 Tutukova S , Tarabykin V , Hernandez-Miranda LR . The role of neurod genes in brain development, function, and disease. Frontiers in Molecular Neuroscience. 2021;14 :13. doi: 10.3389/fnmol.2021.662774 PubMed PMID: WOS:000664304600001. 34177462 76 Crespi EJ , Denver RJ . Roles of stress hormones in food intake regulation in anuran amphibians throughout the life cycle. Comparative Biochemistry & Physiology A-Comparative Physiology. 2005;141 (4 ):381–90. PubMed PMID: ISI:000232093100004. doi: 10.1016/j.cbpb.2004.12.007 16140236 77 Piprek RP , Kubiak JZ . Development of gonads, sex determination, and sex reversal in Xenopus. Kloc M , Kubiak JZ , editors. Malden: Wiley-Blackwell; 2014. 199–214 p. 78 Kanamori A , Brown DD . The regulation of thyroid hormone receptor beta genes by thyroid hormone in Xenopus laevis. J Biol Chem. 1992;267 (2 ):739–45. PubMed PMID: .1730665 79 Tata JR . Autoinduction of nuclear hormone receptors during metamorphosis and its significance. Insect Biochem Mol Biol. 2000;30 (8–9 ):645–51. PubMed PMID: ISI:000088815500003. doi: 10.1016/s0965-1748(00)00035-7 10876107 80 Lemkine GF , Raji A , Alfama G , Turque N , Hassani Z , Alegria-Prevot O , et al . Adult neural stem cell cycling in vivo requires thyroid hormone and its alpha receptor. FASEB Journal. 2005;19 (7 ):863–5. doi: 10.1096/fj.04-2916fje|ISSN0892-6638. PubMed PMID: WOS:000227591900010.15728663 81 Billon N , Tokumoto Y , Forrest D , Raff M . Role of thyroid hormone receptors in timing oligodendrocyte differentiation. Dev Biol. 2001;235 (1 ):110–20. PubMed PMID: ISI:000169701100009. doi: 10.1006/dbio.2001.0293 11412031 82 Morte B , Manzano J , Scanlan TS , Vennstrom B , Bernal J . Aberrant maturation of astrocytes in thyroid hormone receptor alpha 1 knockout mice reveals an interplay between thyroid hormone receptor isoforms. Endocrinology. 2004;145 (3 ):1386–91. doi: 10.1210/en.2003-1123 PubMed PMID: WOS:000189035500049. 14630717 83 Perez-Juste G , Garcia-Silva S , Aranda A . An element in the region responsible for premature termination of transcription mediates repression of c-myc gene expression by thyroid hormone in neuroblastoma cells. J Biol Chem. 2000;275 (2 ):1307–14. PubMed PMID: ISI:000084836600082. doi: 10.1074/jbc.275.2.1307 10625678 84 Bedo G , Pascual A , Aranda A . Early thyroid hormone-induced gene expression changes in N2a-beta neuroblastoma cells. J Mol Neurosci. 2011;45 (2 ):76–86. doi: 10.1007/s12031-010-9389-y PubMed PMID: WOS:000295173600002. 20506002 85 Denver RJ , Ouellet L , Furling D , Kobayashi A , Fujii-Kuriyama Y , Puymirat J . Basic transcription element-binding protein (BTEB) is a thyroid hormone-regulated gene in the developing central nervous system. Evidence for a role in neurite outgrowth. J Biol Chem. 1999;274 (33 ):23128–34. doi: 10.1074/jbc.274.33.23128 .10438482 86 Lebel JM , Dussault JH , Puymirat J . Overexpression of the beta-1 thyroid receptor induces differentiation in Neuro-2a cells. Proc Natl Acad Sci U S A. 1994;91 (7 ):2644–8. PubMed PMID: ISI:A1994ND60200053. doi: 10.1073/pnas.91.7.2644 8146169 87 Gudernatsch JF . Feeding experiments on tadpoles. I. The influence of specific organs given as food on growth and differentiation. A contribution to the knowledge of organs with internal secretion. Wilhelm Roux Arch Entwicklungsmech Org. 1912;35 :457–83. 88 Leloup J , Buscaglia M . Triiodothyronine, hormone of amphibian metamorphosis. Comptes Rendus Hebdomadaires Des Seances De L Academie Des Sciences Serie D. 1977;284 (22 ):2261–3. PubMed PMID: ISI:A1977DL25700018. 89 Krain LP , Denver RJ . Developmental expression and hormonal regulation of glucocorticoid and thyroid hormone receptors during metamorphosis in Xenopus laevis. Journal of Endocrinology. 2004;181 (1 ):91–104. PubMed PMID: ISI:000221056200009. doi: 10.1677/joe.0.1810091 15072570 90 Hu F , Knoedler JR , Denver RJ . A mechanism to enhance cellular responsivity to hormone action: Kruppel-like factor 9 promotes thyroid hormone receptor-beta autoinduction during postembryonic brain development. Endocrinology. 2016;157 (4 ):1683–93. doi: 10.1210/en.2015-1980 .26886257 91 Sun GH , Fu LZ , Shi YB . Epigenetic regulation of thyroid hormone-induced adult intestinal stem cell development during anuran metamorphosis. Cell and Bioscience. 2014;4 :8. doi: 10.1186/2045-3701-4-73 PubMed PMID: WOS:000348468000001. 24507416 92 Fu LZ , Yin J , Shi YB . Involvement of epigenetic modifications in thyroid hormone-dependent formation of adult intestinal stem cells during amphibian metamorphosis. General and Comparative Endocrinology. 2019;271 :91–6. doi: 10.1016/j.ygcen.2018.11.012 PubMed PMID: WOS:000456353800009. 30472386 93 Tuinhof R , Ubink R , Tanaka S , Atzori C , van Strien FJ , Roubos EW . Distribution of pro-opiomelanocortin and its peptide end products in the brain and hypophysis of the aquatic toad, Xenopus laevis. Cell Tissue Res. 1998;292 (2 ):251–65. doi: 10.1007/s004410051056 .9560468 94 Yao M , Westphal N , Denver R . Distribution and acute stressor-induced activation of corticotrophin-releasing hormone neurones in the central nervous system of Xenopus laevis. Journal of Neuroendocrinology. 2004;16 (11 ):880–93. doi: 10.1111/j.1365-2826.2004.01246.x. PubMed PMID: WOS:000226255500002. 15584929