Genomic evidence for the parallel regression of melatonin synthesis and signaling pathways in placental mammals

Background: The study of regressive evolution has yielded a wealth of examples where the underlying genes bear molecular signatures of trait degradation, such as pseudogenization or deletion. Typically, it appears that such disrupted genes are limited to the function of the regressed trait, whereas pleiotropic genes tend to be maintained by natural selection to support their myriad purposes. One such set of pleiotropic genes is involved in the synthesis ( AANAT, ASMT) and signaling ( MTNR1A, MTNR1B) of melatonin, a hormone secreted by the vertebrate pineal gland. Melatonin provides a signal of environmental darkness, thereby influencing the circadian and circannual rhythmicity of numerous physiological traits. Therefore, the complete loss of a pineal gland and the underlying melatonin pathway genes seems likely to be maladaptive, unless compensated by extrapineal sources of melatonin. Methods: We examined AANAT, ASMT, MTNR1A and MTNR1B in 123 vertebrate species, including pineal-less placental mammals and crocodylians. We searched for inactivating mutations and modelled selective pressures (dN/dS) to test whether the genes remain functionally intact. Results: We report that crocodylians retain intact melatonin genes and express AANAT and ASMT in their eyes, whereas all four genes have been repeatedly inactivated in the pineal-less xenarthrans, pangolins, sirenians, and whales. Furthermore, colugos have lost these genes, and several lineages of subterranean mammals have partial melatonin pathway dysfunction. These results are supported by the presence of shared inactivating mutations across clades and analyses of selection pressure based on the ratio of non-synonymous to synonymous substitutions (dN/dS), suggesting extended periods of relaxed selection on these genes. Conclusions: The losses of melatonin synthesis and signaling date to tens of millions of years ago in several lineages of placental mammals, raising questions about the evolutionary resilience of pleiotropic genes, and the causes and consequences of losing melatonin pathways in these species.


Plain language summary
Evolution is typically thought to occur by making an organism more complex, such as through the addition of new traits. However, evolution can also proceed through the loss or degeneration of characteristics. The reduction of eyes in animals living underground and the loss of limbs in snakes are examples of such regressive evolution. When organisms evolve regressively, the genes underlying such traits often become broken (pseudogenes) or deleted from their genome altogether. However, genes are not typically thought to be lost when they perform many functions (pleiotropy). For instance, four genes are involved in the production (AANAT, ASMT) and detection (MTNR1A, MTNR1B) of melatonin, a hormone that signals an animal's organ systems that it is dark outside (e.g., during the night). Melatonin genes seem unlikely to become broken or lost given that this hormone is responsible for signalling many parts of the body to perform functions specific to nighttime. The loss of such genes would likely have a negative, cascading effect on many body functions; therefore, natural selection would presumably retain functional versions of these genes. Surprisingly, however, certain vertebrates are reported to lack the organ responsible for secreting melatonin (pineal gland). We examined the melatonin genes in these vertebrates to determine if they are broken, missing and/or show evidence of degradation over time. We discovered that while pineal-less crocodylians retain all four melatonin genes, other species reported to lack a pineal gland (anteaters, sloths, armadillos, pangolins, dugong, manatee, whales) have broken and/or missing melatonin genes. Furthermore, colugos (flying lemurs) and some mammals that live underground in complete darkness show evidence of melatonin gene dysfunction. Together, these results indicate that these mammals lost their ability to produce and detect melatonin tens of millions of years ago, raising questions as to how they have adapted to a life without this hormone.

Introduction
Evidence for the molecular basis of regressive evolution, or vestigialization, has become abundant following increases in the availability of whole genome assemblies ( The pineal gland is an endocrine organ located within the diencephalon of vertebrates, whose primary, and perhaps only, function is to secrete melatonin. Melatonin functions as a potent antioxidant and can also act as a hormone that signals environmental darkness (Hardeland & Poeggeler, 2003;Zhao et al., 2019). In vertebrates, melatonin derives from serotonin, which is converted into N-Acetylserotonin by aralkylamine N-acetyltransferase (AANAT) and then modified into melatonin by N-Acetylserotonin O-methyltransferase (ASMT/HIOMT). Melatonin synthesis in the pineal gland is under the control of the circadian master clock (suprachiasmatic nucleus), following a pattern of high production in darkness and low production in light. It is released into the bloodstream and arrives at target tissues to activate downstream pathways via the G protein-coupled melatonin receptors type 1A (MTNR1A/ MT1) and 1B (MTNR1B/MT2) (Cipolla-Neto et al., 2014). These melatonin receptors are expressed widely, including in the suprachiasmatic nucleus, thalamus, cerebral cortex, retina, kidneys, adrenal glands, reproductive organs, arteries, immune cells, liver, pancreas, skin and bone, indicating broad signaling from this hormone (Cipolla-Neto & Do Amaral, 2018;

Amendments from Version 1
This manuscript has been revised in a few important ways: 1. we elaborated on the methodology provided for the evolutionary analyses (specifically dN/dS analyses), 2. we performed new analyses on historical evolutionary pressures of crocodylian melatonin genes, 3. we provided additional supporting information validating the pseudogene status of nearly every putatively inactivated melatonin pathway gene in our focal taxa. Beyond this, we modified text and added supplementary information to provide further clarity to our study. None of the results changed in any significant way from what was reported in the first published manuscript. If anything, our conclusions have been strengthened by the addition of new data and analyses.
Any further responses from the reviewers can be found at the end of the article

REVISED
Sapède & Cau, 2013). Indeed, by providing a signal of darkness, melatonin helps to maintain circadian and circannual rhythms, influencing energy metabolism, seasonal reproduction, migration behavior, blood pressure, immune system functioning, among other processes ( Given the ubiquity and myriad effects of this hormone, we hypothesized that melatonin synthesis has been maintained even in pineal-less vertebrates, but may rely on extra-pineal sources to perform the same functions. We set out to determine if the genes necessary for melatonin synthesis (AANAT, ASMT) and signaling (MTNR1A, MTNR1B) are functionally intact in several clades of apparently pineal-less vertebrates. We found evidence that in contrast to crocodylians, which maintain intact melatonin pathway genes, numerous placental mammal lineages show evidence of melatonin synthesis and/or signaling disruption. We inferred that several of these events occurred tens of millions of years ago, raising questions about the evolutionary resilience of pleiotropic systems.  3.5 (Edgar, 2004). We then used the human and chicken assembly-derived sequences as the references for obtaining orthologs in mammals and crocodylians, respectively.

Materials and methods
For whole genome assemblies, we BLASTed our reference sequences against assemblies in external and imported databases (Underlying data, Supplementary Table S1 [Emerling et al., 2021]) using intermediate sensitivity settings (e.g., discontiguous megablast on NCBI). When obtaining AANAT sequences, the whole reference sequence was used as a probe in BLAST searches and short read mapping. ASMT sequences were obtained using a mixture of single exons plus flanking DNA and whole gene reference sequences due to the relatively large size of the gene and small exons. In the case of some ASMT sequences, we queried NCBI's annotated scaffolds directly (RefSeq). MTNR1A and MTNR1B both have two coding exons separated by a large intron, so to avoid incompatible homology issues when attempting to align the introns, we typically designed separate probes with flanking sequences for the two coding exons. If we ever failed to recover sequences for a species, we used mRNA sequences and annotated gene predictions, plus newly-assembled sequences, as reference sequences for BLAST and mapping, especially from close relatives. Sequences derived from genome assemblies that contain long stretches of Ns can cause issues with alignments, so we trimmed any such instances to ten Ns.
Short reads were obtained from both published (NCBI, Sequence Read Archive) and novel sources. To generate novel sequences for multiple crocodylians (Caiman latirostris, Crocodylus niloticus, Mecistops sp., Melanosuchus niger, Osteolaemus tetraspis, Paleosuchus palpebrosus, P. trigonatus, Tomistoma schlegelii), high quality DNA was extracted from tissue using the Qiagen DNeasy kit and quantified through Qubit. In preparation for library construction, 200ng of DNA were sheared with a Covaris S220 for an average size of 500 bp and used as input for the Illumina Neoprep automated library construction instrument. Constructed libraries were pooled and sent to the New York Genome Center for sequencing on HiSeq X system paired-end 150 bp reads.
We imported short reads into Geneious and mapped them to reference sequences (exon + flanking) that were taxonomically close (i.e., same genus, family or order) using intermediate sensitivity settings (medium-low sensitivity). For each set of probes, we performed an initial mapping run to expedite sequence capture, followed by a remapping of the captured reads with the "Fine Tuning" option set to iterate up to five times. Mapped short reads were examined by eye and trimmed of nonhomologous sequences likely derived from adapters or sequencing errors.
For dN/dS analyses in crocodylians, we obtained sauropsid outgroup comparison sequences derived from predicted gene models only (Underlying data, Supplementary Tables S2-S5 [Emerling et al., 2021]), not constructing any such outgroup sequences from assemblies. We BLASTed (discontiguous megablast) each Alligator mississippiensis ortholog against NCBI's nucleotide collection, restricting searches to Squamata, Testudinata and Aves, respectively.
Upon obtaining the mammalian and sauropsid gene sequences, each sequence was imported into Geneious and aligned to its respective probe (i.e., human or chicken whole gene sequence) using MUSCLE. When working with different clades, alignments for each species in that clade were performed successively to provide anchoring and improve subsequent alignments, followed by manual adjustment (see clades in the underlying data, Supplementary  Finally, following the completion of our analyses, a whole genome assembly (WGS) was released for the sirenian Steller's sea cow (Hydrodamalis gigas). The high similarity to other sirenian gene sequences (pairwise comparison with Dugong dugon: AANAT: 98%; ASMT: 97.4%; MTNR1A: 96%; MTNR1B: 96.6%) allowed us to positively identify these genes and characterize their functionality, but we did not include them in the evolutionary analyses.

Evolutionary analyses
Following examination of each clade-specific alignment, global alignments for each gene were assembled to perform phylogenetic (RAxML) and selection pressure (PAML) analyses, respectively (Underlying data, Supplementary  In some cases, relationships for certain mammals were not resolved in the reference trees, but confamilials and/or congeners were present, allowing us to confidently place such taxa in the phylogeny. For example, our primary reference tree does not include Desmodus rotundus, but other phyllostomid species are present. Therefore, D. rotundus was positioned in the phyllostomid portion of the PAML tree topology. After setting the topology, we executed one-ratio branch models with one of three codon frequency models (1/61 each, F1X4, F3X4) and used the Akaike information criterion to select the best-fitting codon frequency model for each gene (Underlying data, Supplementary To test whether inferred pseudogenes are under relaxed selection, we employed branch tests on a series of nested models relative to a 'master' model to estimate whether ω is elevated on relevant branches. To construct the 'master' models, we used an approach similar to that employed by Meredith et al. (2009), categorizing branches as either functional, pseudogenic, mixed/transitional or pre-mutation. 'Functional' branches are those for which there is no evidence for gene inactivation, and there is no expectation that the gene is likely to become a pseudogene. All 'functional' branches were grouped together for a single ω, were expected to be consistent with purifying selection (ω < 1) and therefore were treated as the background with which to compare the remaining three branch categories. 'Pseudogenic' branches are those which post-date a branch on which a gene is inferred to have become a pseudogene. For example, if we mapped an inactivating mutation to the stem cetacean branch, all of the descendant branches within crown Cetacea were grouped together, and a single ω was estimated for this set of branches. In such cases, the dN/dS is expected to approach or be under relaxed selection (ω = 1, higher than functional category). 'Mixed/transitional' branches are those on which a gene is inferred to have become a pseudogene, due to one or more inactivating mutations being mapped to that branch. It is therefore suspected to have a mixed history, transitioning from a functional gene to a pseudogene. Depending on how early or late this transition has occurred, the ω is expected to be highly similar to the functional category (more recent inactivation), very similar to the pseudogene category (very early inactivation) or in between these two extremes. Finally, the 'pre-mutation' branch category represents instances in which we found no positive evidence of pseudogenization on a branch, but external evidence suggests that relaxed selection may have occurred. Such pre-mutation branches could be designated due to multiple immediate descendant lineages having inactivating mutations and/or one or more other key genes being inactivated. As an example of the former, if a gene is inferred to have been pseudogenized within sloths, armadillos and anteaters separately, given that these all belong to the xenarthran clade, we would designate the branches prior to the mixed/transitional branches of these subclades as 'pre-mutation' branches. As an example of the latter, if a taxon had one or both melatonin synthesis genes (e.g., Fukomys damarensis) and/or both melatonin receptor genes (e.g., Heterocephalus glaber) inactivated, but no evidence of pseudogenization in the remaining genes, we would categorize these branches as 'pre-mutation'.
As we were not interested in average selection patterns across all branches within a particular category, except the background (functional) branches, we allowed each pre-mutation branch, each mixed/transitional branch and each set of pseudogene branches that post-dated an inactivation event to have their own ω estimates. This allowed us to test whether unique historical instances of putative gene inactivation were consistent with relaxed selection, rather than obscuring the signal by 'averaging' pre-mutation, mixed/transitional and pseudogene branches, respectively, across the tree.
After setting these branch categories, this 'master' model (Extended data, Supplementary Figures S1-S4 [Emerling et al., 2021]), was calculated for each gene and used as the basis for which to compare subsequent nested models (Underlying data, Supplementary Tables S7-S10 [Emerling et al., 2021])). Specifically, if higher than the background in the master model, we tested whether each set of pseudogene category branches, each mixed/transitional branch and each pre-mutation branch was distinguishable from the background (functional) ω and/or from the neutral ω value of 1 (expectation for pseudogene branches). In such cases, we used the master model for each gene that included all branch categories (Extended data, Supplementary Figures S1-S4 [Emerling et al., 2021]), and would change only one branch of interest to be fixed as being part of the background or 1. We then compared each of the nested models to the master model using likelihood ratio tests to determine which models better fit the data for a branch of interest. In these instances, the models in which branch(es) were fixed as the background or 1 are considered the null model, with the master model being the alternative model, i.e., we tested whether a free estimate of ω for a foreground branch is statistically distinguishable from the background and/or 1. Finally

Melatonin gene expression
In order to test for evidence of active transcription of the melatonin synthesis genes in a crocodylian, we analyzed published RNA sequencing data from 22 sequencing experiments on tissues in a juvenile Alligator mississippiensis (St John et al., 2012). We BLASTed (megablast) the protein-coding regions of the A. mississippiensis AANAT and ASMT genes against short reads derived from mRNA sequencing in NCBI's Sequence Read Archive (Underlying data, Supplementary Table S12 [Emerling et al., 2021]), imported the sequences into Geneious, and mapped them to the A. mississippiensis references using the low sensitivity setting.

Results
We found evidence that all xenarthrans, all pangolins, and all cetaceans have had both their melatonin synthesis (AANAT, ASMT) and melatonin signaling (MTNR1A, MTNR1B) genes disrupted via the accumulation of inactivating mutations and whole gene deletions (  "Melatonin synthesis disrupted" indicates AANAT and/or ASMT is inferred to have been inactivated on the associated branch. "Melatonin signaling disrupted" indicates MTNR1A and MTNR1B are both inferred to have been inactivated on the associated branch, or one gene was lost on an earlier branch and the second was inactivated on the associated branch. Note that the stars are arbitrarily placed in the middle of branches and do not correspond to a precise timing for gene loss. Letters on stars and nodes correspond to letters in

Xenarthra
For AANAT, sloths (Folivora) share a mutated start codon (ACA), splice donor mutation (GT to GG) in intron 1, 1-bp deletion in exon 2, and two frameshift indels in exon 3. All anteaters (Vermilingua) possess a start codon mutation (GTG or TTG), with six additional inactivating mutations shared between Myrmecophaga and Tamandua. Armadillos (Cingulata) share a 1-bp deletion in exon 1, and four frameshift indels in exon 3 ( Figure 2). In ASMT, there is a putative 1-bp insertion in exon 6 shared by two-toed sloths (Choloepus spp.), two chlamyphorid armadillos (Calyptophractus, Chlamyphorus), and a dasypodid armadillo (Dasypus novemcinctus French Guiana). This region is missing in all other xenarthran sequences (i.e., deleted or not assembled), and the alignment is admittedly ambiguous. Furthermore, in all xenarthrans we examined, exon 5 appears to be absent and the ancestral stop codon (i.e., the stop codon common to all other species examined) is mutated, although it differs among sloths (TAT), anteaters (TGT or TGG), and armadillos (CGT or TGT). Beyond these, shared inactivating mutations are present among the anteaters (5-bp deletion, exon 2 [ Figure 2]; two frameshift indels, exon 8), sloths (1-bp deletion, exon 3; 4-bp insertion, exon 7), and armadillos (1-bp deletion, exon 6), respectively. We were unable to obtain MTNR1A in most xenarthrans, possibly due to a whole gene deletion, although the assemblies were not complete enough to verify this via a synteny analysis. Nonetheless, the silky anteater (Cyclopes didactylus) and sloths retain one or more MTNR1A pseudogenes (Figure 2). Sloths appear to have multiple paralogs, the identities of which are difficult to tease out, with separate exons being found on separate contigs, but the pattern suggests all are probably inactivated. For MTNR1B, a 2-bp deletion in exon 2 is shared by sloths and anteaters (Pilosa; Figure 2). Among armadillos, a single 1-bp deletion (exon 1) is shared by dasypodids, and five inactivating mutations are shared among chlamyphorid armadillos (two in exon 1; three in exon 2).
From dN/dS ratio analyses, we found evidence of shifts in selection pressures consistent with ancient pseudogenization of these genes. Estimates in which key xenarthran branches had statistically elevated ω values were found for AANAT ( . We also ran a model for ASMT in which crown Xenarthra was given a single ω, under the assumption that it was inactivated in the stem lineage (see above). This ω was estimated as 0.79, and when compared to a model in which the crown Xenarthra ω was fixed at 1, the former was a better fit for the data than the latter. As such, this result seems to be inconsistent with a stem xenarthran inactivation of ASMT.

Pholidota
For pangolins, all three species share two premature stop codons in AANAT (exons 1 and 2), another in exon 5 of ASMT, and four inactivating mutations in MTNR1B (Figure 2). Exon 1 of MTNR1A is missing for both Manis spp., preventing comparisons with Phataginus tricuspis, but the former two share a 1-bp deletion in exon 2 ( Figure 2)

Dermoptera
Both melatonin synthesis genes have shared inactivating mutations across all three colugos: exon 1 of AANAT and exon 3 of ASMT each have shared 1-bp deletions, (Figure 2). We were unable to obtain MTNR1A for Cynocephalus volans, but both Galeopterus species have a splice donor mutation (AG to AT) in the intron (Figure 2). Similarly, we were unable to assemble exon 2 of MTNR1B for C. volans, and exon 1 appears intact, but Galeopterus spp. share a 1-bp insertion (exon 1; Figure 2 Other placental mammals Beyond these taxa, AANAT is present as two to three paralogs in multiple non-cetacean cetartiodactyls, with one paralog sometimes being a pseudogene, a finding corroborated by a recent study (Yin et al., 2021). However, at least one gene is always intact in all non-cetacean cetartiodactyls that we examined. ASMT is a pseudogene in three subterranean rodents (Fukomys damarensis [Bathyergidae]; Nannospalax galili, Spalax carmeli [Spalacidae]), MTNR1A is inactivated in two hyraxes (Procavia capensis, Heterohyrax brucei, Hyracoidea), a monk seal (Neomonachus schauinslandi), pig (Sus scrofa), talpid mole (Condylura cristata), and the naked-mole rat (Heterocephalus glaber), and MTNR1B is a pseudogene or likely deleted (i.e., negative BLAST results) in a host of other species, including seven afrotherians, eight carnivorans, one bat, all six examined eulipotyphlans, one lagomorph, five primates, and five rodents we examined ( Figure 1). The few examples of shared inactivating mutations we found in these other species include a 4-bp insertion in exon 6 of ASMT in two spalacids (Extended data, Supplementary Figure S9  Here also, dN/dS ratio analyses suggest that many of these pseudogenes are indeed under relaxed selection based on statistically elevated ω values.

Melatonin genes inactivated in many mammals, but functional in crocodylians
Here we reported evidence that xenarthrans, pangolins, cetaceans and some sirenians have lost the capability to synthesize and bind melatonin via the traditional pathway found in vertebrates, coinciding with the ostensible absence of a pineal gland in these taxa. This builds upon recent studies on cetaceans and the West Indian manatee (Huelsmann et al., 2019; Lopes-Marques et al., 2019), demonstrating the surprising extent of the degradation of these genes. We hypothesized that despite the pineal gland's apparent absence in xenarthrans, pangolins, sirenians, cetaceans and crocodylians, the genes underlying the production and signaling of melatonin would remain intact, given the widespread effects of melatonin in vertebrates and evidence of circulating melatonin in these taxa  , 1980), this suggests that either the pineal gland is intact but difficult to isolate and/or extra-pineal sources of melatonin exist. In support of the latter hypothesis, we found evidence of AANAT and ASMT expression in the eyes of the American alligator (Alligator mississippiensis). Furthermore, a previous study found that the gene encoding a pineal opsin pigment is a pseudogene in crocodylians (Emerling, 2017b), potentially revealing a shift in the source of melatonin from the pineal gland to the eye, following the regression of the former.  1985). A recent anatomical study failed to find a pineal gland in the white-bellied pangolin (Phataginus tricuspis), although the absence was attributed to an error in preparing the brain for study (Imam et al., 2019). Despite these inconsistencies, a regressed, if not completely absent, pineal gland largely predicts the disappearance of the canonical melatonin pathways, at least in mammals. Furthermore, we found evidence of melatonin synthesis inactivation in colugos, suggesting that these species may also possess a regressed pineal gland, or may lack it entirely. Although one study found a diencephalon in colugos comparable in relative size to gliding rodents and bats (Pirlot & Kamiya, 1982), a review of pinealocytes in mammals noted the lack of research specifically on the pineal in Dermoptera (Bhatnagar, 1992).  (Hardeland & Poeggeler, 2003). However, the functional significance of serum melatonin may be obviated by the absence of melatonin receptors in these species, presumably preventing any contribution to circadian signaling in the body's tissues. Despite this, given that recent discoveries have shown that presenting as a pseudogene does not always indicate that all biological function is lost (Cheetham et al., 2020), there remains a possibility that these apparently dysfunctional genes are able to contribute to melatonin metabolism in some unknown fashion.

Ancient loss of melatonin genes in placental mammals
Our results strongly suggest that melatonin synthesis and signaling have been abolished within multiple placental mammal lineages for extensive periods of evolutionary time. The most ancient of these may be in Xenarthra, given our evidence of possible ASMT and MTNR1A inactivation on the stem branch. Crown Xenarthra arose roughly 68 million years ago (mya) (Gibb et al., 2016), potentially meaning that the loss of melatonin synthesis and possibly some signaling took place near the K/Pg boundary, when non-avian dinosaurs went extinct and placental mammals began to radiate. However, given some ambiguity in the evidence for stem inactivation of ASMT (see results), and the absence of sequences of MTNR1A for most species of xenarthrans, convergent loss remains a strong possibility. Despite this, additional shared mutations in AANAT and ASMT suggest that the components for melatonin synthesis were disrupted prior to the origin of armadillos (45 mya), sloths (31 mya), and anteaters (38 mya), respectively. In addition, MTNR1B was likely pseudogenized prior to the sloth / anteater split (59 mya), and the origin of chlamyphorid armadillos  et al., 2019), in order to help modulate circadian and circannual physiological processes, melatonin synthesis and signaling would seemingly become indispensable for most vertebrates. As such, it is a challenge to clarify the causes and consequences of losing melatonin pathway genes (Valente et al., 2021). Convergent evolution often results from similar selection pressures, which may explain why both cetaceans and some sirenians have lost these genes. Perhaps the unique demands of a fully aquatic lifestyle, such as the need to frequently surface for respiration, needed to be uncoupled from a rhythmic signal of darkness. By contrast, the semi-aquatic sea otter (Enhydra lutris) and pinnipeds have retained their melatonin synthesis genes, although MTNR1B is pseudogenized in E. lutris and two seals (Phocidae), with one of the phocids also showing evidence of MTNR1A inactivation. Perhaps this underlies their intermediate aquatic phenotype, although eight of the 13 carnivorans also present an MTNR1B pseudogene. Another example of strong convergent evolution can be seen in xenarthran anteaters and pangolins, both of which have radically modified their feeding apparatus to ingest ants and termites (myrmecophagy). However, other myrmecophagous taxa we examined, including the aardvark (Orycteropus afer), aardwolf (Proteles cristatus), and bat-eared fox (Otocyon megalotis) at most only have MTNR1B inactivated.
Regardless of their specific phenotypes, all of these taxa experience fluctuations in light and darkness, so it is unclear as to why loss of such a hormone would be beneficial. By contrast, it seems more logical for melatonin synthesis to be lost while adapting to an environment of nearly complete darkness, in which rhythmic secretions entrained on light patterns may no longer be possible. Multiple lineages of subterranean mammals fit this description, and indeed, we found evidence of ASMT pseudogenes in the subterranean rodents Nannospalax galili, Spalax carmeli (Spalacidae) and Fukomys damarensis (Bathyergidae), and inactivation of both receptors in the naked mole-rat (Heterocephalus glaber; Bathyergidae). The latter result had been previously reported in a single individual (Kim et al., 2011), but we confirmed that both genes share the same disabling mutations in a second individual H. glaber and are likely under relaxed selection. The fossorial star-nosed mole, Condylura cristata (Talpidae), also appears to have dispensed of both melatonin receptors, with MTNR1A a pseudogene and MTNR1B being completely absent from the assembly. In addition, dN/dS estimates suggest selection on ASMT is relaxed in H. glaber, a species which also has an atrophied pineal gland (Quay, 1981), and MTNR1B is inactivated in N. galili, F. damarensis and a golden mole (Chrysochloris asiatica; Chrysochloridae). This may be relevant to the ancient loss of melatonin synthesis and signaling in xenarthrans, given that comparative anatomy and an analysis of genes critical for vision in bright light appear to point to an early subterranean history for this clade (Emerling & Springer, 2015). Perhaps an extended history underground limited the utility of melatonin synthesis and signaling, and upon emerging from this committed existence in the darkness, their descendants inherited this unusual phenotype. However, a pineal gland that can synthesize melatonin has been reported in at least one spalacid (Balemans et al., 1980), pineal glands are reported to be present in talpids and chrysochlorids (Legait et al., 1976;Pevet, 1974;Pevet & Kuyper, 1978) and the melatonin synthesis genes remain intact for Chrysochloris asiatica and Condylura cristata.
One potentially unifying concept for the pattern of pineal gland / melatonin synthesis loss may be related to thermoregulation. Ralph (1975) hypothesized that the size of the pineal gland in vertebrates may be correlated to thermoregulation, tentatively linking pineal gland size to latitude, activity pattern and relative homeothermy. He observed that some of the largest pineal glands belong to species that inhabit higher latitudes, while pointing out that some vertebrates with the smallest or absent pineal glands tend to be restricted to the tropics. Though the comparisons were limited, a better-substantiated pattern was noted with the related parietal eye of squamates. The parietal eye is an organ that is developmentally related to and anatomically linked with the pineal gland, which appears to largely provide information for thermoregulation in ectothermic squamates, possibly through melatonin regulation (Firth & Kennaway, 1980;Hutchison & Kosh, 1974;Phillips & Harlow, 1981;Stebbins & Eakin, 1958). When comparing the presence or absence of the parietal eye, researchers noted a trend of parietal eye loss in squamates that live near the equator (Gundy et al., 1975). Given that the parietal eye provides information about temperature and light, which largely correlates with the amount of sunlight, and the pineal gland secretes melatonin in darkness, the latitudinal hypothesis may have some validity. At lower latitudes, there is less seasonality; therefore, being able to detect changes in circadian and circannual dark and light cues is plausibly of less adaptive benefit in these regions. Notably, xenarthrans, pangolins, sirenians and colugos live almost exclusively in the tropics; furthermore, aquatic and subterranean habitats provide a buffering effect from temperature fluctuation. These characteristics encompass all taxa we record as lacking melatonin synthesis capabilities, making it a potentially attractive hypothesis.
Significantly, the patterns of melatonin pathway degradation have strong overlap with placental mammals that have lost the capacity for non-shivering thermogenesis (NST). Specifically, xenarthrans, pangolins, cetaceans and sirenians all have inactivated UCP1, a gene that facilitates NST (Gaudry et al.,  2017). Furthermore, hyraxes (Hyracoidea) and pigs (Suidae) have a UCP1 pseudogene, hyraxes (Procavia capensis, Heterohyrax brucei) have both melatonin receptor genes inactivated, and the wild boar (Sus scrofa) has an MTNR1A pseudogene. Notably, melatonin induces the production of brown adipose tissue (Heldmaier et al., 1981; Heldmaier & Hoffmann, 1974), a major location of NST. Together, these data suggest that the loss of melatonin synthesis may be coupled with the loss of this thermoregulatory tool in these clades, further underscoring a potential link between melatonin pathway loss and changes in thermoregulatory requirements.

Conclusion
In this study we provided evidence that, in contrast to crocodylians, numerous placental mammals reported to lack a pineal gland also lack the genes necessary for the canonical vertebrate melatonin synthesis and signaling pathways. However, this result seems to raise more questions than answers. Given the pleiotropic nature of melatonin synthesis and signaling genes, which selection pressure(s) could have led to the loss of this seemingly crucial signaling molecule? What are the physiological consequences of this loss? Are there possibly compensatory alternative mechanisms for producing and sensing melatonin? For those species that present serum melatonin, how are they doing so? Does this melatonin function merely as an antioxidant, or does it aid in circadian and circannual signaling via different pathways? Further studies on comparative anatomy, physiology and gene expression in pineal-less taxa and others should shed further light on these challenging questions.

Reviewer Report 14 December 2021
Many other lineages also showed evidence of inactivation or loss in a smaller number of the genes. The authors speculate that loss of the melatonin synthesis and receptor genes may be associated with changes in thermoregulation and the loss of non-shivering thermogenesis. Overall, the manuscript is very well written, presents new and interesting results, and outlines several important areas for future research. I outline below several areas where the authors could further improve the paper: Inference of gene loss and pseudogenization depends on high quality sequencing data and genomic assemblies and apparent loss and pseudogenization could be due to incomplete sequencing/assembly or sequencing/assembly errors. While patterns of shared indels or premature stop codons are particularly convincing, I believe this caveat should still be addressed, especially considering that ~70% of the mammals examined were listed as having at least one pseudogenized or lost gene.

○
Further to the point above, the sources of the genomic data used in the study is unclear.
The main supplementary table of the genome sources does not provide much meaningful information (e.g., source listed as WGS). Looking through some of the subsequent SI tables accession numbers are given for some of the genes, but others are simply listed as New.
Ideally all of the genomes would be made available, but minimally more detail needs to be provided on the genomic data (quality, completeness, etc…) used to obtain the sequences analyzed in the current study. It is difficult to evaluate the robustness of the results or reproduce the study without this.

○
It was also unclear whether the transcriptome data were newly sequenced or the analysis was of previously published data. More detail should be provided.

○
The explanation of the methodology and approach used in the selection analyses was difficult to follow and perhaps overly complex. The authors mention a nested approach but do not provide much detail. Looking through the supplementary tables of PAML results I do not see nested sets of models or an attempt to find the best fitting models. From what I can tell based on the tables, the approach relies on what appears to be a free ratio model, the use of which is generally discouraged due to it being overly parameter rich. A large number of likelihood ratio tests are then made between the free ratio model and models placing individual branches or clades in the 'foreground'. This approach appears to (incorrectly) use the most parameter-rich model as the null model for the likelihood ratio tests. Models that restrict the 'foreground' branch/clade to ω = 1 are also used but they are also compared to the parameter rich model (rather than to same model but where the foreground ω is free to vary). As far as I can tell, no models are used that combine multiple branches or clades as might be expected. Beyond these issues, a much simpler approach could be used that also more directly tests the author's hypotheses by comparing ω between mammals with and without pineal glands and between lineages with and without evidence for pseudogenization in that specific gene. A codon model specifically designed to detect relaxation of selective constraint (ie, RELAX) would be useful in this regard.
○ It would be nice to see selection analyses on the crocodilian sequences as well to, for example, show that they are still under strong selective constraint. I think this could still be useful without other sequences for comparison (eg, from birds/turtles). ○ I initially found the mix of analyses on mammalian and crocodilian genes unusual, but this ○ made more sense within the context of animals with and without pineal glands. This could perhaps be made more clear, especially in the title which only mentions placental mammals and not crocodylians or the lack of pineal glands.
I generally thought the figures were well put together quite effective, but the text size (most notably in Fig. 2) can be quite small and hard to read, even when zooming in.

Is the work clearly and accurately presented and does it cite the current literature? Yes
Is the study design appropriate and does the work have academic merit? Yes

If applicable, is the statistical analysis and its interpretation appropriate? Partly
Are all the source data underlying the results available to ensure full reproducibility? Partly We agree that this can indeed be an issue with genomic data. While we are confident that the shared inactivating mutations, often paired with dN/dS estimates suggestive of relaxed selection, are strong evidence of pseudogenization in the pineal-less mammals and colugos, we believe Dr. Schott is correct that we should present additional evidence for the inactivation of genes in other taxa. To address this, we have done the following: Collected additional genomic sequences from related taxa for species with seemingly unique inactivating mutations (i.e., to test for evidence of shared mutations).

1.
In cases where this did not provide positive evidence of inactivating mutations (e.g., no adequate sister taxon, no shared mutations), we mapped short reads to the genomic contigs for such taxa.

2.
We performed these analyses for all taxa that have evidence for inactivation of one or more melatonin synthesis gene and/or both melatonin receptor genes, given that these are the most relevant for our study. We thank Dr. Schott for this critical feedback, as Dr. Janiak had the same concerns; it clearly needed to be clarified further. In short, we did not use a free ratio model as the null model. We designated a model in which certain branches were grouped based on their association with distinct inferences of pineal gland loss / gene losses, with branches being defined as functional (no evidence of dysfunction), pseudogenic (postdating branch in which pseudogenization is inferred to have occurred), mixed/transitional (branch in which pseudogenization is inferred to have occurred), and pre-mutation (no evidence of pseudogenization on the branch but external evidence suggests gene may be under relaxed selection). As stated in the manuscript, this is following the approach of Meredith et al. 2009. Accordingly, there was a 'master model' for each gene, with each branch category described above, and this model was then compared to null models in which a specific pseudogene branch (or set of branches), mixed/transitional branch or pre-mutation branch was set to the background or 1, effectively testing whether a specific branch or set of branches was distinct from the background and/or 1, and therefore displaying evidence of relaxed selection. This was to test for individual historical cases of relaxed selection, rather than a broad 'average' pattern of relaxed selection, the latter of which does not accomplish what we were trying to test. Details of what we added can be found in the sections described below: We have elaborated on the textual specifics of these methods under a section entitled "Evolutionary analyses" within the Materials and Methods (pgs. 9-10). Figures S1-S4). We then paired the branch labeling with the dN/dS results in Supplementary Tables S7-S10 by use of numbering, to make for a smoother comparison of the two.

2.
We believe that this should be far easier to understand than in our previous iteration, and hope it will satisfy Drs. Schott and Janiak.
It would be nice to see selection analyses on the crocodilian sequences as well to, for example, show that they are still under strong selective constraint. I think this could still be useful without other sequences for comparison (eg, from birds/turtles).
We agree with this suggestion and have done so. We added sections in the Materials and Methods (pgs 8 & 10) and Results (pg 13), added the alignments and topology to Supplementary Datasets S5-S9, and modelling results in Supplementary Tables S6 and S11. The basic result is that all four genes appear to be under purifying selection in crocodylians, with three ratio models (background, stem Crocodylia, crown Crocodylia) being statistically indistinguishable from one ratio models, and the crown and stem Crocodylia w estimates being very similar to the background w.
I initially found the mix of analyses on mammalian and crocodilian genes unusual, but this made more sense within the context of animals with and without pineal glands. This could perhaps be made more clear, especially in the title which only mentions placental mammals and not crocodylians or the lack of pineal glands.
Regarding the title, given that the first version of the manuscript is technically published via the Open Research Europe model, we thought changing it would make it confusing for future referencing purposes, especially as the manuscript has already been cited. Furthermore, given that the crocodylians we examined did not show evidence of inactivating melatonin pathway genes, including a crocodylian reference in the title seemed to us to take away from our most important finding, i.e., mentioning the absence of positive evidence for inactivation in crocs may detract from the title. However, we have tried to make the inclusion of crocodylians in our study clearer in the following ways: Including an initial mention of the crocodylian results at the end of the introduction (pg 6).
secretes melatonin. Given the importance of melatonin in regulating an organism's circadian rhythm, among other function, the authors hypothesized that the four genes involved in melatonin synthesis and signaling may remain functional even in species that lack a pineal gland, as there may be other sources of melatonin within the body.
Finally, I am admittedly not an expert in melatonin pathway or circadian rhythm research, so I hope that others can provide more insightful feedback on the authors' discussion of the implications of their findings. In my non-expert opinion the authors appear appropriately cautious in how they interpret their findings, posing more questions than answers.

3.
A statement in the Materials and Methods on the generation of crocodylian sequences (on pg. 7). 1.
An explicit statement that the RNA sequencing data for the Alligator has been previously published with a reference to the corresponding paper (pg. 10).

2.
The letters along the branches in Figure 1 are quite small and hard to read without zooming in. Could these be a bit bigger?
We have addressed this by making the letters larger and superimposed colored circles upon to make them stand out among the branches.
Finally, I am admittedly not an expert in melatonin pathway or circadian rhythm research, so I hope that others can provide more insightful feedback on the authors' discussion of the implications of their findings. In my non-expert opinion the authors appear appropriately cautious in how they interpret their findings, posing more questions than answers.
We agree that this would have been ideal. We requested reviews from three different melatonin experts, but unfortunately none were available. Due to time constraints, we have decided to revise the manuscript with the reviews provided so far and will look forward to commentary from future studies.

Competing Interests:
No competing interests were disclosed.