Thyroid hormones (THs) are essential for the normal growth, development and health of all vertebrates (1–3) Among many other biological processes, THs control metabolism, cell differentiation and migration and also play a pivotal role for e.g., neurodevelopment (3, 4). One of the most remarkable examples of TH action is anuran metamorphosis, where a tadpole transitions into a frog with extensive reprogramming and reorganization of organs and tissues (5–8). This remodeling has direct parallels to the developing human during the last trimester and the first weeks of life, where a surge of THs is essential for proper development (8–10). Hence, anuran metamorphosis is an excellent model system to study developmental endocrinology and the dynamics of the TH system. This has been acknowledged in several publications (5, 7–9, 11). The potential of amphibian metamorphosis as a screening tool for chemical disruption of the hypothalamus-pituitary-thyroid axis has furthermore been recognized in the OECD guideline 231 for the testing of chemicals i.e., the Amphibian Metamorphosis Assay (12).

Despite this, the current knowledge of baseline levels of THs during metamorphosis largely relies on pooled samples using methodologies known for high variability with sparse reporting of measured variation (13–23) ( Table 1 ). In fact, (radio)immunoassays are still one of the most commonly used methods for quantification of THs despite known limitations, e.g., cross-reactivity and interference from proteins as summarized by Spencer (24). The current study uses a highly robust, selective and sensitive liquid chromatography tandem mass spectrometry (LC-MS/MS) method based on isotopic dilution to quantify 12 THs and metabolites (16, 25).

THs exert their chemical messenger effect mainly through nuclear TH receptors, although growing evidence shows that non-genomic signaling also plays a role (8, 26). In the genomic pathways, nuclear TH receptors (THRs) bind directly to DNA or indirectly by interacting with chromatin associated proteins to modulate expression levels of hormone-responsive genes (26). Two THR genes, α and β, code for these receptors. Both have been found to show specific expression patterns according to cell type and developmental stage (27). In addition to TH concentrations, the transcription levels of THR genes, therefore, serve as robust indicators of TH signaling and the endocrinological processes during metamorphosis. Earlier studies have shown that an abundance of transcripts encoding THR α and β (thra, thrb) as well as TH-induced bZip protein (thibz) show high correlation with developmental stage in anuran tadpoles (27–29) indicating their potential to serve as bioindicators of metamorphosis in addition to THs. However, these data and their relationships to changes in TH levels in the same subjects have not been established in amphibians.

American bullfrog (Rana [Lithobates] catesbeiana) is a great alternative model species as it is a representative of the largest amphibian group, Ranids or “true frogs”, and is found worldwide. Their large-sized tadpoles enable the analysis of individual animals rather than pooled samples. Furthermore, this amphibian compared to e.g., Xenopus laevis tadpoles, also has a greater similarity to human life stages as they completely transition from aquatic to terrestrial living (Xenopus remain fully aquatic as adults). By using wild caught animals instead of laboratory manipulated individuals, we ensure that the natural variation of endocrinological state is represented. By using a high sample replication, we also capture the natural biological variation without losing statistical power.

Although anuran metamorphosis is initiated solely by THs, many other biochemical pathways undergo extensive remodelling during this transition (30) and the regulation of metamorphosis is a complex process of interacting mechanisms (11). The complexity, tight regulation, and interconnection of biological networks and pathways can be conceptualized using untargeted metabolomics. In contrast to traditional hypothesis-driven investigations, untargeted approaches, such as metabolomics, interrogate the whole biological system and let the organism-response guide the investigations (31–35). Such knowledge is very powerful when trying to make new discoveries about a biological system and can aid the search for new robust biological markers. In the case of anuran metamorphosis, this type of analysis provides a holistic view of the biological processes and allows for a search of robust bioindicators of the molecular dynamics during the transitions. Such bioindicators are useful for assessing the state of the endocrine system and its disruption by external environmental factors (such as endocrine disrupting chemicals or climate change (i.e., temperature)).

The first aim of the present study is to map the dynamics of known bioindicators, ie. TH concentration and transcription of TH sensitive genes in wild-caught R. catesbeiana tadpoles and frogs during metamorphosis. This includes the quantification of measurable THs and metabolites in serum as well as abundance of gene transcripts of thra, thrb, and thibz in back skin of the same animals. As a second aim, we use untargeted metabolomics to reveal new bioindicators of metamorphosis.

The present study brings forth essential new knowledge of baseline levels of THs and TH-responsive genes during metamorphosis. Currently, published data largely rely on pooled samples, limited replication, and sparse reporting of variation ( Table 1 ). In contrast, the present study includes individual sampling of 5-10 animals per developmental stage and the methodology reflects state-of-the-art quantification of biomolecules. This also includes the quantification of TH metabolites other than T3 and T4, such as T2, which provides novel insight into the thyroid hormone system. Table 1 summarizes maximal levels of THs (ng/mL) measured during anuran metamorphosis in published studies. The average reported value for all species is 6.7 ± 0.6 ng/mL and 2.0 ± 0.6 ng/mL for T4 and T3, respectively. The studies of R. catesbeiana show slightly lower average values of 6.0 ± 0.9 ng/mL and 1.1 ± 0.9 ng/mL for T4 and T3, respectively. In contrast, the concentrations measured in the present study are higher. As shown in the final row of Table 1 , we report maximum levels of 28 ± 10 ng/mL and 2.3 ± 0.5 ng/mL for T4 and T3, respectively. The reason for this dissimilarity could reflect differences in the developmental stage that was sampled for quantification (see Gosner stage in Table 1 ) and the masking of variability by using pooled samples. Additionally, the extraction method can cause differences especially for T4. T4 is largely protein-bound (51), and to obtain a quantification of total T4, the extraction method in the present study included a protein denaturation by urea to release all protein-bound hormone. This differs from earlier applied methodologies and could explain the higher quantity measured in the present study. Notably in the only two previous studies examining T4 levels in individual tadpoles, it was noted that all stages had some individuals where T4 was not detected (20, 21). In the present study, we found Gosner stage 44 animals all had detectable T4 concentrations but two individuals had T4 levels below the LoQ. The present study results are, therefore, consistent with these earlier findings in that T4 is highly variable, which contributes to the high variation in T4 levels during metamorphic climax. The average fold change during metamorphosis is furthermore much greater for T4 i.e. 20 fold, compared to the 5 fold change in T3, which explains why measurements in T4 spans a concentration range all the way from below the detection limit at Gosner stage 31-39 up to the maximum measurement of 64 ng/mL at Gosner stage 44.

In addition to updated information on circulating thyroid hormones and metabolites (T4, T3, T2, T0 and T1Am), this is also the first study to simultaneously analyse same-individual back skin gene expression of the thyroid hormone receptors α and β (thra, thrb) and TH-induced bZip protein (thibz) throughout metamorphosis. Efforts have been made earlier to combine published data on X. laevis to get this simultaneous view (8, 52), but the same individual measurement, with high level of replication offers an unprecedented depiction of the endocrinological dynamics during anuran postembryonic development, and by approximation, vertebrate development as well. In previously published illustrations, the peak in T3 concentration was reported to occur simultaneously with the peak in thrb transcript abundance and precede the maximum in T4 (23, 52). In the present study, however, we observed that the peak in T3 and T4 occurs simultaneously, and that thrb abundance preceded the increased concentration of both T3 and T4. This exemplifies the exquisite timing of events during the metamorphic process as the synthesis of the receptor occurs in anticipation of the peak in TH concentration.

The observed timewise correlation between THs and the abundance of thibz and thrb was confirmed by correlation analyses. Furthermore, the PCA showed that T3 and T4 as well as thrb and thibz transcript abundance were the main factors to separate the metamorphic climax from other developmental stages. Because the back skin transcript data are so well linked to the TH measurements it opens up the possibility of using non-lethal skin swabs on the rostrum area as a validated proxy for TH levels to detect the status and the possible disruption of the TH system. The current results are however species- and tissue specific, so this would, at this point, still have to be further verified on a case-by-case basis.

The abundance of thra transcripts did not, however, follow the same pattern; its abundance remained stable until the end of metamorphic climax where it decreased significantly and remained low in the juveniles. This significant reduction coincides with maximal T4/T3 and it could potentially be an example of negative feedback. The observed metamorphosis-related dynamics is indeed in good correspondence with an earlier study in X. laevis, where a plateau in thra expression also was observed until a drop after metamorphic climax (53). Relative THR expression is, however, known to be tissue specific, so it is important to note that our results illustrate the transcript abundance in back skin, a tissue that remodels but is retained throughout metamorphosis. Other tissues may exhibit different dynamics as exemplified by the differences in tissue responsiveness to THs (29) and tissue specific expression profiles of THRs (54). In particular, thra seems to have a very diverse expression pattern across tissues (27, 29, 54).

To the best of our knowledge, this is the first study to quantify T2 throughout tadpole metamorphosis. Interestingly, T2 increased earlier than T3 and T4 prior to the significant increase in transcript abundance of thrb and thibz and remained high during metamorphosis. This indicates a specific and distinct biological role of this biomolecule in anuran postembryonic development. There is limited research on the action of T2 in different vertebrates and, to the best of our knowledge, there are no studies about the biological role in amphibians. In humans and rodents it has, however, been established that this iodinated molecule has biological activity distinct from that of T4 and T3 (55), including for instance activation of fatty acid oxidation and down-regulation of lipogenesis in hepatic fatty acid metabolism (56–61). In the freshwater fish, tilapia, T2 acts as a selective activator of a specific long isoform of THR-b named L-TRb1 and modulates its expression in vivo (58). T2 also specifically regulates gene networks associated with cell signaling and transcriptional pathways in brain and liver, which are distinct from those affected by T3 (59, 60). Furthermore, T2 affects myelination in cerebellum (59), which is essential for neurodevelopment and participates in growth processes via L-TRb1 activation (61).

Based on these findings and the fact that T2 concentration was significantly increased earlier than T3 and T4, it is tempting to speculate that T2 may modulate THR expression in amphibians. More research is clearly needed to investigate this hypothesis, including the expression and activity of deiodinases. Generally, TH-related research tends to focus on T4 and T3, of which the biological role is well-documented. However, deiodinated products like T1 and T2, as well as aminated metabolites (e.g., T1Am) and iodothyroacetic acids (e.g., 5-diiodothyroacetic acid, Diac), are increasingly being recognized as biologically active (57, 62, 63). These molecules will be an important focus in furthering our understanding of TH dynamics.

The present work on bullfrog tadpole postembryonic development may facilitate the understanding of TH-dependent processes during the perinatal period of human development. The role of THs in the transition from an “aquatic” environment to a “terrestrial” one affects multiple organ systems that are paralleled in human and bullfrog (8–10). It is interesting to speculate that T2 in addition to T3 and T4 plays a more prominent role and remains to be investigated.

The untargeted analysis showed that the changes in THs during metamorphosis are accompanied by similar dynamics in many other biomolecules. Earlier untargeted studies in the serum of R. catesbeiana using ultraperformance liquid chromatography coupled with mass spectrometry (UPLC-MS) have found a high degree of remodelling in five core metabolic pathways during metamorphosis. This included arginine and purine/pyrimidine, cysteine/methionine, sphingolipid, and eicosanoid metabolism and the urea cycle (30). In addition, a major role for lipids during the postembryonic process was established and arachidonic acid metabolism also underwent remodeling during metamorphosis (30). In accordance with this, the current study also observed metamorphosis-related changes in several lipids including eicosanoids and arachidonic acid metabolites ( Figure 3 ). The bioindicators that showed the highest predictive power for distinguishing metabolite dynamics like the THs, however, included prostaglandins, acyl-carnitines and fatty acyls ( Table 2 ).

L-Carnitine is an endogenous metabolite found in all mammalian species (64). The most important role of L-carnitine is to facilitate the transport of fatty acids into the mitochondria for beta-oxidation (64). During this process, L-carnitine forms a high-energy bond with a fatty acid to give an acyl carnitine derivative that is transported into the mitochondrial matrix. Here it undergoes transesterification with coenzyme A (CoA) to form the corresponding acyl-CoA, which then enters the fatty acid b-oxidation (64). In this way, L-carnitine is also able to regulate the intracellular ratio of free CoA to acyl-CoA by forming acylcarnitines and eventually remove surplus of acyl groups from the organism by preferential renal excretion. This means that the intramitochondrial relationship between free CoA and acyl-CoA is reflected in the extracellular L-carnitine to acylcarnitine ratio, giving indirect evidence of altered mitochondrial metabolism (64). Considering this, our observation of altered carnitine levels during metamorphosis, which additionally correlates with the TH-dynamics, is consistent with the fact that the mitochondrial metabolism is modified during metamorphosis. A tissue specific metabolic reorganization during metamorphosis has already been observed for R. omeimontis tadpoles (65). This study found that, by metamorphic climax, these tadpoles showed overall increased energy metabolism. In the liver, this was accompanied by upregulated amino acid catabolism and downregulated β-oxidation and glycolysis, suggesting a metabolic shift from carbohydrate and lipids to amino acid. A similar shift is likely occurring in R. catesbeiana. Therefore, the results in the current study show that the acyl-carnitines are excellent bioindicators for the metabolic shift during metamorphosis. Further cross-species validation is however needed.

The ratio of L-carnitine to acyl-carnitines is currently applied in neonatal screening for inborn errors of metabolism, so the potency of these bioindicators is already known from other vertebrates (66). Acyl carnitine concentrations and THs can be biochemically linked as T3 is a known inducer of carnitine palmitoyltransferase I gene expression in the liver (67). This enzyme catalyzes the transfer of long chain fatty acyl groups to carnitine from CoA for translocation across the inner mitochondrial membrane. A long chain fatty acyl, which was tentatively annotated as 2R-HpOTrE, and its oxidation product, 9-OxoOTrE, were actually also among the 20 compounds selected on sPLS-DA component 1. The effect of THs on the transferases for short- and medium chained fatty acyls to the best of our knowledge, not been investigated. However, the two metabolites that were classified as medium-chain hydroxy acids and derivatives ( Table 2 ) could potentially be related to this.

We also observed modifications to the arachidonic acid metabolism pathway and eicosanoids during metamorphosis ( Figure 3 ). Furthermore, a prostaglandin was among the 20 compounds that best discriminate between Gosner stages ( Figure 3C , Table 2 ). Prostaglandins are synthesized by sequential enzymatic modification of arachidonic acid (68). They are mainly recognized for their role in inflammatory responses (69) but they also play a role in neurogenesis (70). A link to TH homeostasis have been observed in rats with high iodide intake-induced hypothyroidism where increased levels of prostaglandins were found (71). Therefore, prostaglandins and other eicosanoids present an additional chemical class of important bioindicators of metamorphosis and vertebrate developmental endocrinology that could be applied in future studies.

The steroid hormone, testosterone, was interestingly found to follow a similar Gosner stage profile as T3, T4, carnitines and prostaglandins. Testosterone is a steroid hormone and a Cyp17-derivative of pregnenolone or progesterone precursors. Cyp17 expression have recently been found in gonad-mesonephros complex in X. laevis tadpoles, showing that androgen synthesis is possible (72). Testosterone is the main sex hormone in vertebrate males. However, because the animals in the study were of mixed sex, testosterone seems to serve a different purpose in this developmental window. Although we cannot exclude tissue-specific transformation of testosterone to estrogens through aromatase activity, it should be noted that the observed increase in testosterone occurred after anticipated sexual differentiation and is completely independent of sexual maturation. Much empirical evidence exists that confirms crosstalk between the hypothalamic-pituitary-thyroid (HPT) axis which involves TH synthesis and the hypothalamic-pituitary-gonadal axis which involves synthesis of steroid hormone (7, 73). As reviewed in Thambirajah et al. (7), this includes the effect of androgens like testosterone on the HPT axis. The exact mechanisms and physiological consequences are still unclear; however, androgens have been linked to bone formation and growth (74), which may be the reason for the observed metabolic pattern during metamorphosis.

The vital crosstalk between the HPT axis and the hypothalamic-pituitary-adrenal axis, which signals the production of corticosteroids, is also well established. This includes hormones like cortisol and corticosterone. Corticosterone is maximal at metamorphic climax much like T3 and T4 and it can accelerate TH-dependent metamorphosis. It has also recently been shown to be required for survival through metamorphosis in X. tropicalis (75). A metabolite tentatively annotated as corticosterone was found among the compounds that displayed same Gosner stage profile as T3 and T4, however, it could not be annotated with a spectral match. These findings underline what many others has called for before, namely that crosstalk between endocrine axes receives more attention in future studies. It is one more piece in the puzzle that is essential for our understanding of vertebrate developmental endocrinology. The non-target data provides a necessary context for interpreting TH-mediated developmental processes and potential molecular crosstalk mechanisms. Furthermore, this provides a solid foundation from which alternative bioindicators of metamorphosis may be chosen as feasible endpoints in future risk assessment of chemicals and their potency for causing developmental toxicity.

The datasets presented in this study can be found in online repositories. The names of the repository and accession number can be found below: https://massive.ucsd.edu/ProteoSAFe/static/massive.jsp, MSV000090714.

Treatment of animals was in accordance with the guidelines of the Canadian Council on Animal Care and experimental protocol #2015-028 was approved by the University of Victoria Animal Care Committee.

RP: Conceptualization, methodology, investigation and formal analysis, writing – original draft. SJ: Methodology, investigation and formal analysis, writing – review and editing. MH: Conceptualization, methodology, writing – review and editing, funding acquisition and resources. CH: Conceptualization, methodology, investigation and formal analysis, writing – review and editing, funding acquisition and resources. All authors contributed to the article and approved the submitted version.