Molecular evolution of TRPC4 regulatory sequences supports a role in mammalian thermoregulatory adaptation

Introduction

The transient receptor potential (Trp) ion channels are proteins that form nonselective, gated transmembrane channels, permitting the flow of Ca²⁺ and other ions upon activation by regulating molecules (Ramsey, Delling & Clapham, 2006). In vertebrates, the ‘canonical’ TRP subfamily (TRPC) of proteins participate in diverse processes including cardiac and neurological function. They also function in the transduction of environmental stimuli such as thermoregulation, proprioception, chemoreception and the perception of pain (Minke & Cook, 2002). Major regulatory molecules that are known to alter TRPC channel activity include calmodulin, inositol triphosphate receptors (ITPRs), phospholiplase C, stromal interaction molecule 1 (STIM1), and G-protein coupled receptors, as well as several scaffolding proteins and small molecules such as diacylgycerol and phosphatidylinositol 4,5-bisphosphate (PIP₂). These molecules primarily interact with various C-terminal residues of TRPCs that project interior to the plasma membrane (Kim et al., 2019; Chen et al., 2020).

Trpc4 has been shown by Zhou et al. (2023) to be essential in mice for non-shivering thermogenesis (NST), a suite of behavioral and metabolic responses to core temperature change that are regulated by warm-sensing neurons (WSNs) in the preoptic area of the hypothalamus (Tan et al., 2016). A key mechanism effecting NST is metabolic heat generation by brown adipose tissue (BAT), although thermogenesis can also occur in muscle (Nowack et al., 2017). BAT thermogenesis is mediated by thermogenin, also known as uncoupling protein 1 (UCP1), and is considered an ancestral trait of eutherian mammals (Keipert et al., 2024). BAT thermogenesis is important during acute cold stress (particularly in small-bodied mammals) and for postnatal survival, and may play a role in the regulation of hibernation in some species (Hashimoto et al., 2002; Heldmaier, Klaus & Weisinger, 1990; Bienboire-Frosini et al., 2023).

Trpc4 has separately been identified as a positive-selection candidate in bats (order Chiroptera) and exhibits a higher rate of diversification of the regulatory protein-binding region in that clade compared to related mammalian orders (Cornman, 2024). In light of the functional analyses of Zhou et al. (2023), this increased diversification of Trpc4 in bats may be associated with thermoregulatory divergence within that clade. Numerous studies have documented ecologically important differences in the thermoregulatory physiology of bats, including the use of hibernation and torpor as well as tolerance of changes in core body temperature (Tb) while active (e.g., Cryan & Wolfe 2003; Stawski, Willis & Geiser, 2014; Czenze & Dunbar, 2017; Czenze et al., 2022; Rubalcaba et al., 2022). More generally, NST exhibits strong evolutionary lability in eutherian mammals, as evidenced by the independent pseudogenizations of Ucp1 in several eutherian groups, which has been associated with Pliocene climate change and increased body size (Gaudry et al., 2017). Thus, ecological diversification of mammalian thermoregulation might be reflected in protein evolution of TRPC4 at various phylogenetic scales if such changes are important modulators of NST.

Support for the hypothesis that diversifying selection on TRPC4 regulatory sequences is related to thermoregulatory adaptation in mammals can be evaluated by comparing evolutionary rates of Trpc4 and interacting genes in clades with different thermoregulatory strategies. One possible comparison is among vertebrate lineages prior to and since the advent of NST, for example by comparing predominantly ectothermic taxa (amphibians and nonavian sauropsids), taxa that thermoregulate without BAT-mediated NST (birds), and taxa that thermoregulate using BAT-mediated NST, among other mechanisms (eutherian mammals). This comparison would reveal whether Trpc4 evolutionary rates changed in a manner consistent with functional divergence and subsequent purifying selection due to recruitment into NST-regulating genetic pathways. A complementary comparison is made possible by the independent pseudogenizations of Ucp1 noted above. Assuming the Trpc4-dependence of NST in mice is characteristic of mammals generally, evolutionary rates within the regulatory region of TRPC4 may increase following Ucp1 loss due to the elimination of that functional constraint. An additional comparison is suggested by the pronounced effect of body size on core body temperature variation (Tb) within a species (Boyles et al., 2013). Small mammals experience a wider range of Tb than larger mammals in general, and useful generalizations can be made at the order level as to relative body size and importance of NST without neglecting the importance of adaptations within particular lineages. For example, taxa exhibiting high Tb variation are common within rodents (order Rodentia) and bats, whereas various carnivores (order Carnivora) are prominent among taxa with low variation in Tb (Boyles et al., 2013). These order-level differences allow potential associations between the evolutionary rate of Trpc4 and the level of variation in Tb (heterothermy) to be explored. Importantly, while functional overlap may exist between genetic regulation of NST and the prolonged shifts in Tb caused by hibernation and torpor, there is as yet no basis for predicting a correlation between use of these latter strategies and Trpc4 evolutionary rate.

A complication of this analytical framework is that TRPC channels assemble as tetramers and may be heteromeric (i.e., containing more than one type of TRPC protein) and thus Trpc genes may not evolve independently. TRPC4-containing channels in mouse neurons are predominantly heteromers with TRPC5 or TRPC1 (Bröker-Lai et al., 2017; Kollewe et al., 2022), and Trpc5 expression has itself been linked to physiological responses to temperature (Zimmerman et al., 2011; Bernal et al., 2021). In fact, functionally significant amino-acid variation affecting non-NST temperature response has been identified in the C-terminal regulatory region of Trpc5 (Ptakova et al., 2022), the phylogenetically closest paralog of Trpc4 (Chen et al., 2020). In contrast to Trpc1, Trpc4, and Trpc5, proteins encoded by Trpc gene family members Trpc3, Trpc6, and Trpc7 form predominantly homomeric channels in neurons (Kollewe et al., 2022). Including other Trpc genes in a comparative analysis would therefore allow patterns unique to Trpc4 to be differentiated from patterns that are common to many genes compared across a given phylogeny. The same reasoning applies to structurally similar genes of other Trp subfamilies that have identified roles in transducing temperature signals in various contexts. Trpm2 is of particular note as it acts as a homeostatic regulator of elevated temperature (as opposed to reduced temperature) in warm-sensitive neurons (Kamm et al., 2021; Song et al., 2016). Recent work has implicated adiposal expression of Trpa1, Trpv2, and Trpm2 in the control of BAT thermogenesis at the source (Sun, Uchida & Tominaga, 2017; Lv et al., 2023; Benzi et al., 2024), but roles in the homeostatic activation or maintenance of NST responses within warm-sensitive neurons themselves have not been reported. Thus, while the present study hypothesizes neuronal control of evolutionary variation in NST phenotypes, manifested by changes in TRPC4 in heteromeric channels, regulatory changes in thermogenic tissues themselves are a reasonable alternative.

In summary, this study investigates Trpc4 evolutionary rates with respect to NST at different phylogenetic scales. Evolutionary rates are estimated for a discrete, unambiguously defined regulatory region of sufficient length and conservation to be informative of selection pressures acting on it. Where statistical support of rate heterogeneity is found, additional statistical models are applied to differentiate positive selection and relaxed constraint as explanatory factors. Trpc4 evolutionary rates are then interpreted within a larger context of interacting and non-interacting Trpc homologs as well as structurally similar Trp family members with potentially overlapping functions.

Methods

Predicted transcripts were downloaded from the ortholog browser of the National Center for Biotechnology Information (NCBI) for each gene and each taxonomic group. NCBI ortholog assignments are performed for high-quality reference genomes based on reciprocal sequence similarity to model genes as well as microsyntney of local gene order (refer to https://www.ncbi.nlm.nih.gov/kis/info/how-are-orthologs-calculated/), and are considered authoritative for the objectives of this study. The major phylogenetic bins within vertebrates for which sequence was accessed were amphibians, avian sauropsids (class Aves), non-avian sauropsids (lizards, turtles, snakes, and their relatives), eutherian mammals, and marsupials. Orthologs for two monotreme species (Ornithorhynchus anatinus and Tachyglossus aculeatus) were also accessed but after inspection were not combined with marsupials for most analyses because the two groups aligned poorly for several genes, which would have substantially reduced the number of codons available for analysis. However, monotremes were included in the branch-specific test of Trpc4 divergence in basal Eutheria (see Results), as they aligned unambiguously and represented important variation for rate estimation. Within eutherian mammals, sequences were downloaded for the orders Chiroptera (bats), Artiodactyla (even-toed ungulates and cetaceans, synonymous with Cetartiodactyla (sensu Hassanin et al., 2012), Primates, Rodentia (rodents), and Carnivora (carnivores). These are the major mammalian orders for which the number of high-quality gene annotations are highest (refer to Genes database of NCBI), and thus best-suited for order-level evolutionary rate comparisons. Fish were not included in the analysis as most Trpc genes are duplicated within that clade, consistent with ancestral polyploidization (Von Niederhäusern et al., 2013). Due to fluctuation in taxonomic nomenclature used by NCBI, the following species pairs are considered synonymous: Neogale vison and Neovison vison (synonymous also with Mustela vison); Lagenorhynchus obliquidens and Sagmatias obliquidens; Microtus fortis and Alexandromys fortis; Physeter catodon and Physeter macrocephalus.

Nucleotide alignments were initially generated with MAFFT v. 7.490 (Katoh et al., 2002) and edited with BioEdit v. 7.7.1 (Hall, 1999) by first trimming upstream of start codons using the coordinates annotated by NCBI for each transcript. Alignments were then collapsed, translated, and re-aligned at the amino-acid level with ClustalW v. 1.81 (Thompson, Gibson & Higgins, 2003) before reverting to nucleotide sequence to trim alignments downstream of stop codons. Ambiguously aligned regions (e.g., gapped sequence, low-complexity sequence, or splice-site variants) were manually trimmed between the nearest unambiguously aligned codons, and sequences that were substantially incomplete in the studied region were removed. Analyzed alignments are available in File S1, and by policy are also made available in a USGS data release (Cornman, 2025).

Gene family members were dispersed in the studied genomes (i.e., they do not occur as tandem arrays that may be more prone to mis-annotation or mis-assembly) and no indication of gene conversion between paralogs was observed in the alignments (i.e., tracts of sequence similar between paralogs in the same taxon but dissimilar to orthologs from related taxa). All Trp gene alignments were confirmed to begin with the conserved TRP domain (the “TRP Box” depicted in Fig. 1.4 of Rosasco & Gordon, 2017) that marks the transition between transmembrane channel domains and the intracellular tail of the protein (Duan et al., 2018; Vinayagam et al., 2018). The initial residues of the domain are highly conserved within each Trp family (Rosasco & Gordon, 2017) and readily identifiable. Regulatory binding sites, alternative splice sites, and other sequence features of the region were extracted from (Tang et al., 2001; Mery et al., 2002; Odell, Van Helden & Scott, 2008; Miehe et al., 2010; Yuan et al., 2003) as well as the Genome Viewer utility of NCBI (Rangwala et al., 2021). As previously reported by Frankenberg et al. (2011), Trpc2 was found to be inconsistently annotated due to its close genomic proximity to and past conflation with Xndc1 and was not analyzed here. A full-length alignment and dendrogram of these proteins were made for selected species to further confirm their orthology (Fig. S1). The sequences were amino-acid aligned with Mega v. 11 (Kumar et al., 2018) and a neighbor-joining tree was generated using the JTT exchange matrix and the default gamma distribution of rate variation.

Draft guide trees for evolutionary rate analysis were generated using approximately 15 kb of mitochondrial sequence downloaded from NCBI for available taxa. Mitochondrial sequences were aligned without regard to gene features with MAFFT and then trimmed of the poorly aligned control region. Trees were generated using the neighbor-joining algorithm in Mega with the maximum composite likelihood substitution model and with rate variation among sites modeled by the default gamma distribution. Taxa used in Trp gene alignments that lacked sufficient mitochondrial sequence were inserted manually into Newick-formatted trees based on established taxonomic classifications (Hassanin et al., 2012; Agnarsson et al., 2011; Blanga-Kanfi et al., 2009; Hassanin et al., 2021; Zurano et al., 2019). Guide trees were unrooted with ape (Paradis, Claude & Strimmer, 2004) prior to use and are available in File S1.

Protein evolutionary rates ω were estimated with the codeml program of PAML v. 4.9j (Yang, 2007) under several codon models. To qualitatively compare trends across genes for the various tetrapod clades and mammalian orders, single ω estimates were calculated for each gene-clade combination without rate variation among codon sites or among branches (NSsites = 0, Model = 0). For analyses explicitly parameterizing branch-level rate variation within an alignment, rate classes were marked for each branch in the tree (all branches with the same class designation share a common ω estimate) and model settings changed to NSsites = 0 and Model = 2. The statistical significance of branch-level rate variation for specific comparisons was based on the relative likelihoods of the null model, which had one or two branch rates depending on the hypothesis tested, and an alternative model with one additional rate class (see Results for details). The test statistic was twice the difference in ln(likelihood) of the two models assuming a χ² distribution with one degree of freedom (Yang, 2007), with α set to 0.05. Note that codeml likelihoods are not based on a single most-likely ancestral sequence reconstruction but instead are weighted across all ancestral sequence probabilities reconstructed by the codeml algorithm (Yang, 1998; Anisimova, 2012); ancestral sequences are not specified by the user.

Significant changes in ω under PAML models (see Results) were further evaluated with evolutionary rate models of the HyPhy package (Pond, Frost & Muse, 2004) that explicitly model proportions of sites assigned to rate classes representing positive selection or relaxed selection (i.e., branch-site models). aBSREL (Smith et al., 2015) was used to explicitly test a model of positive selection on Trpc4 in basal Eutheria, as distinct from relaxed selection, as it requires one of three ω classes to be greater than one on the foreground branch for a nonzero proportion of sites. RELAX (Wertheim et al., 2015) was used to explicitly model change in constraint in descendent eutherian mammals relative to evolutionary rates estimated for the basal eutherian branch. The null model under this framework estimates the proportion of sites in a user-designated number of ω classes that may be less than or greater than one, but are constant across the phylogeny. The alternative model estimates a ‘relaxation’ or ‘intensification’ factor k that distinguished two different sets of ω estimates for foreground versus background branches of the phylogeny. For k significantly less than one, the foreground ω classes are closer to one than in background (i.e., both purifying and positive selection are weaker); if k is significantly greater than 1, foreground ω classes are further from one compared with background branches (both purifying and positive selection are stronger). RELAX was used to evaluate whether changes in ω following Ucp1-inactivation were more consistent with relaxation or increased positive selection. Both RELAX and BUSTED (Murrell et al., 2015) were used to distinguish relaxed selection from positive diversifying selection on Trpc genes in bats. BUSTED is similar to RELAX in that three ω classes are assumed by default and are the same for all branches under the null model but differ on foreground branches under the tested model. The highest rate class is capped at one under the null model, whereas it may exceed one in the tested model; positive diversifying selection is inferred if the tested model is significantly more likely and a nonzero proportion of sites have an ω > 1. BUSTED differs from aBSREL in that all foreground branches are treated as a single partition of the data whereas aBSREL tests positive diversifying selection on each branch separately and is therefore less useful for testing pervasive change in a large section of a tree. Note RELAX was run with the default number of rate classes (3) for all analyses except the comparison of descendent eutherian mammals to the basal eutherian branch, which used only two because higher numbers resulted in at least one ω > 1 for the basal eutherian branch, in conflict with the aBSREL results (see Results).

Model outputs are summarized in File S2 and the PAML control file is given in File S1. P-values were adjusted for testing multiple genes or alignments for each question, if applicable, as shown in File S2, using the Benjamin-Hochberg option of the p.adjust function in R (R Core Team, 2018).

Results

Trpc4 evolutionary rate was elevated at the origin of NST

To test whether the evolutionary rate of Trpc4 increased in basal eutherian mammals, a branch-specific test was performed with codeml with three branch classes (Fig. 1). In the null model, marsupials and monotremes were assigned one common evolutionary rate (ω₀) and eutherian mammals a second evolutionary rate (ω₁); two rates were modeled instead of one to reflect the variation among those clades (see below). The alternative (tested) hypothesis specified that the branch leading to eutherian mammals evolved at a third rate (ω₂) that was greater than the other two rates. The data strongly supported an increased divergence rate along the basal branch of Eutheria (P = 4.271E−07), with ω₂ = 0.349 compared to ω₀ = 0.144 and ω₁ = 0.0785. An explicit test of positive selection on this branch (i.e., a nonzero proportion of sites with ω > 1) conducted with aBSREL was not significant (P = 0.210), however. Thus, the data confirm an increased evolutionary rate of TRPC4 regulatory sequences coincident with the evolution of BAT in Eutheria but do not definitively demonstrate the increase was due to positive selection rather than relaxed selection.

Trpc4 evolutionary rates became more constrained after the advent of NST

Evolutionary rates of the Trpc4 C-terminal regulatory region were similar in ectothermic clades and endothermic clades lacking NST, whereas ω was substantially lower in eutherian mammals (Fig. 2). The ω estimates for amphibians and non-avian sauropsids were 0.184 and 0.226, respectively, and 0.187 in birds. In eutherian mammals, ω was estimated to be 0.080, whereas the ω estimate for marsupials was intermediate (0.129). Evolutionary rates of other Trp genes confirm that the above pattern is specific to Trpc4. Only Trpc7 had a lower ω in eutherian mammals than in the other clades, and in that comparison marsupials had the highest ω rather than ectothermic clades. Note Trpc7 was not analyzed for amphibians as there were insufficient orthologous sequences identified for this group.

RELAX found significant evidence of intensification of selection on Trpc4 in descendent eutherian mammals relative to the basal eutherian branch (blue branches relative to green branch in Fig. 1, P = 4.10E−06). A two-rate model estimated an intensification parameter of 2.71 in descendent eutherian mammals, with ω₁ = 0.047 for a proportion of sites p₁ = 0.933 and ω₂ = 1 for a proportion p₂ = 0.067 of sites. These values were ω₁ = 0.324 and ω₂ = 1 in the basal eutherian node (rate proportions are held constant in the model). This explicit statistical test supports the qualitative conclusion drawn from the above comparison of multiple Trp genes across multiple tetrapod groups, i.e., that Trpc4 experienced greater constraint in descendent eutherian mammals after the evolution of BAT in that clade.

Trpc4 evolutionary rates increased after Ucp1 inactivation

Taxa with presumed inactivating Ucp1 mutations include suids, equids, cetaceans, elephantids, manatees, pangolins, armadillos, and sloths (Gaudry et al., 2017). Marking representatives of these taxa as foreground and related mammals as background (Fig. 3) yielded strong support from codeml for increased diversification rates of the TRPC4 regulatory region after Ucp1 inactivation (FDR corrected P = 2.30E−05). The foreground ω was estimated to be 0.101 and the background ω was 0.0574. Note the number of cetacean taxa included in this test was reduced to avoid including very similar sequences, as short branch lengths tend to reduce the power of selection tests (Anisimova, Bielawski & Yang, 2001; Smith et al., 2015).

There are potential concerns with a global test of the effect of Ucp1 inactivation on Trpc4 evolutionary rates, however, as it may not be realistic to assume a single background and a single foreground rate across the entire tree (see Results below). Therefore, the data were split into two additional PAML tests to further confirm the initial result. One test was performed for Boreoeutheria and a second test for “Atlantogenata” sensu Waddell et al. (1999), i.e., Xenarthra plus Afrotheria (Fig. 3). These secondary tests were both statistically significant after FDR correction (P = 0.0116 for Boreoeutheria and P = 0.00210 for Atlantogenata). Within Boreoeutheria, the background ω was estimated as 0.0631 and the foreground ω for Ucp1-inactivated taxa was estimated to be 0.0961. Within Atlantogenata, the background ω estimate was 0.0652 and the foreground ω esimtae for Ucp1-inactivated taxa was 0.119. These parameter values are similar to those estimated for the combined tree.

Relaxed selection on the regulatory region of TRPC4 after Ucp1 inactivation, as opposed to episodic positive selection, was supported by RELAX (FDR-corrected P = 4.49E−06). The relaxed selection model was also significant for the Atlantogenata and Boreoeutheria alignments separately (FDR-corrected P = 5.93E−04 and 6.27E-04, respectively).

Tetrapod TRPC4 sequences have diversified within a distinct, functionally uncharacterized region

While a protein alignment is not informative of evolutionary rate, it can be useful to visualize how extant variation among clades relates to known functional elements within the C-terminus of TRPC4 (Fig. 4). Most of the amino-acid divergence among the tetrapod groups analyzed lies within an ∼85 amino-acid region immediately downstream of a conserved but alternatively spliced sequence (discussed further below) and immediately upstream of the terminal, conserved ‘PDZ’ domain that functions to localize TRPC channels in membranes (Mery et al., 2002). No specific protein-binding interactions have yet been localized to this region of TRPC4. In contrast, most previously characterized protein-binding sites are well conserved among all tetrapod clades examined (Fig. 4). One exception is the second calmodulin binding site identified by Tang et al. (2001), which is less conserved in the ectothermic clades and shows some variability in birds and mammals as well. The region between the second calmodulin binding site and the variable region is believed to bind ITPRs (Tang et al., 2001) but this interaction has not been further localized. These patterns indicate that the variable region upstream of the PDZ binding site is the most likely to be involved in among-clade differences in Trpc4 function, rather than previously characterized binding sites.

Increased evolutionary rate of Trpc4 in bats is not generalizable to small bodied, heterothermic clades

Of the five eutherian orders analyzed, the evolutionary rate of Trpc4 was highest in bats (Fig. 5). Bats also had the highest evolutionary rates for the Trpc4-interacting homologs Trpc5 and Trpc1. Estimates of ω for all three genes were significantly higher in bats than background, represented by outgroups Carnivora and Artiodacytla, in codeml branch tests of evolutionary rate (FDR-adjusted P = 9.64E−08 for Trpc4, P = 2.092E−13 for Trpc5, and P = 0.00361 for Trpc1). Bats did not have higher evolutionary rates for Trpc genes that do not interact with Trpc4 in neurons (Trpc3, Trpc6, and Trpc7), nor did they have elevated rates for the other Trp genes implicated in temperature perception and BAT modulation (Trpm2, Trpv2, and Trpa1). In contrast to bats, ω was estimated to be lower in rodents than in primates for eight of nine genes evaluated, including for Trpc4. Thus, while bats show increased diversification rates of Trpc4 as well as the interacting genes Trpc5 and Trpc1 since their divergence from other mammals, the pattern is not generalizable to other clades exhibiting high variation in Tb.

Analysis with HyPhy favored an interpretation of relaxed selection on Trpc4 and Trpc1 in bats and positive diversifying selection on Trpc5 in bats. RELAX supported relaxed selection on Trpc4 and Trpc1 (FDR-corrected P = 1.590E−11 and 9.13E-4, respectively) whereas the positive selection model of BUSTED was not significant for either. Conversely, the BUSTED model of positive selection within bats was significant for Trpc5 at P = 0.0459, and the RELAX model found evidence of intensified selection on that gene in bats (4.890E-12).

Discussion

NST is a key mammalian innovation but exhibits evolutionary lability as mammalian taxa have adapted to diverse environments and niches (Rodríguez-Serrano & Bozinovic, 2009; Boyles et al., 2013). For example, taxa tolerate different ranges of Tb variation during normal activity and differ in how much Tb can be increased by NST (e.g., Nespolo, Opazo & Bozinovic, 2001; Gębczyński & Taylor, 2004), although species also exhibit substantial plasticity in these traits as well (Van Sant & Hammond, 2008 (Panaino et al., 2023). NST may also be a component of thermoregulatory strategies such as torpor and hibernation, which have distinct endocrinal and metabolic drivers (Hashimoto et al., 2002). Thus, the genetic mechanisms by which NST responses are modulated within and across taxa are likely diverse. Yet the discovery that Trpc4 is required for mouse NST (Zhou et al., 2023) combined with the observation of high diversification rates in the regulatory region of TRPC4 within bats (Cornman, 2024) suggested a hypothesis that relatively simple sequence changes in TRPC4 could underlie aspects of thermoregulatory adaptation in that clade and perhaps in mammals generally.

In this study, multiple comparisons of evolutionary rate were made to assess whether evolutionary change in the regulatory region of TRPC4 was consistent with this hypothesis, which was supported for mammals generally. First, the evolutionary rate of Trpc4 substantially increased after the divergence of eutherian mammals from non-eutherian ancestors coincident with the evolution of BAT (Fig. 1), although no evidence of a rate class ω >1 was found with aBSREL. Trpc4 subsequently became significantly more constrained in descendent Eutheria relative to the basal branch. Moreover, Trpc4 evolutionary rates are substantially lower in eutherian orders than in all clades lacking NST, a pattern not mirrored in any other Trp gene analyzed (Fig. 2), including genes that control Tb homeostasis in warm-sensitive neurons but are not known to function in the induction of NST (Song et al., 2016; Kamm et al., 2021). Additionally, evolutionary rates of Trpc4 significantly increased under relaxed selection in eutherian lineages that later lost thermogenin (UCP1), a major component of NST (Fig. 3). While these results are correlative and no specific amino-acid substitutions in TRPC4 have yet been demonstrated to modify NST regulation by WSNs, the data are nonetheless consistent across multiple independent predictions.

In contrast, no simple relation was identified between heterothermy in small-bodied mammals and Trpc4 evolutionary rates (Fig. 5). Trpc4 as well as Trpc5 and Trpc1 diversified at significantly higher rates in bats than in outgroups, indeed higher than all other eutherian groups analyzed. In contrast, Trpc4 did not have a higher evolutionary rate in rodents compared with primates. Thus, the increased diversification of a suite of Trpc genes that interact in neurons appears to be distinctive of bats and not a recurring pattern in small-bodied, heterothermic eutherian mammals. In fact, Trpc4 and Trpc1 appear to be under relaxed selection rather than diversifying selection in bats, suggesting that regulation of NST by WSNs may be less constrained in species employing frequent daily torpor or prolonged bouts of heterothermy while not actively foraging, as is common in this clade (Stawski, Willis & Geiser, 2014). Alternatively, given the support for positive selection on Trpc5 in bats and the physical interaction of these gene products in heteromeric channels, it is possible that the increased evolutionary rates of Trpc4 and Trpc1 are compensatory to adaptive changes in TRPC5 (e.g., Haag & Molla, 2005; Haag, 2007).

Trpc4, Trpc5, and Trpc1 function in diverse tissues and contexts, often as part of a network of overlapping signals that are integrated at the cellular level (Chen et al., 2020), such that many phenotypic targets of selection could be imagined. Even so, a particularly interesting alternative to thermoregulation is the role of these genes in spatial memory, which is a key aspect of bat behavior (Geva-Sagiv et al., 2015). Triple knockout mice lacking Trpc4, Trpc5, and Trpc1 in hippocampal neurons have deficits in spatial memory and learning (Bröker-Lai et al., 2017). Spatial memory is also impaired in mice that are deficient in TRPC1 alone, either genetically or via pharmacological inhibition (Lepannetier et al., 2018). Furthermore, Trpc4 and Trpc5 are particularly densely expressed in pyramidal neurons of the hippocampus (Fowler et al., 2007), which have long been recognized as key contributors to spatial memory (Burgess, Maguire & O’Keefe, 2002; Geva-Sagiv et al., 2015), and the TRPC4β isoform has been shown to influence the dendritic morphology of hippocampal neurons (Jeon et al., 2013).

Other behavioral phenotypes may also be relevant, as Trpc4 has been linked to autism in multiple human association studies (Gupta et al., 2023; Drago, Calabro & Crisafulli, 2024) and a mouse Trpc4 knockout model shows strong deficits in sociality (Seo et al., 2024). Both Trpc4 and Trpc5 nulls independently diminish anxious or fearful behaviours mediated by the amygdala (Riccio et al., 2009; Riccio et al., 2014), as do pharmacological inhibitors of their protein products (Just et al., 2018). Spontaneous human variants in Trpc5, including partial deletions and nonsynonymous substitutions, are associated with obesity, food hoarding, and impaired sociality in clinical cases and in UK Biobank data (Li et al., 2024). These phenotypes are strikingly recapitulated in mouse models of some of these human variants (Li et al., 2024). Interestingly, two of the nonsynonymous substitutions identified by Li et al. (2024) in the C-terminal regulatory region of human TRPC5 influenced protein decay rates in cell culture, suggesting a possible functional relevance of positive selection on these regions.

While it is increasingly feasible to express modified alleles in cell culture or in genetic model organisms, functional assessment of diverse wild-type alleles in an ecologically relevant context is more challenging, particularly as the molecular pathways that link Trpc4 to the initiation or maintenance of NST are not yet delineated. Nonetheless, combinatorial knock-in studies in mice could provide insights into the ecological significance of sequence divergence in TRPC1/TRPC4/TRPC5 heteromeric channels in bats, as might pharmacological treatment of bats with targeted inhibitors of those proteins (Rubaiy et al., 2017; Just et al., 2018). Comparative assessments of the relative expression of Trpc homologs at the transcript and protein levels (e.g., Li et al., 2022; Lv et al., 2023) could elucidate functional aspects of ecological differentiation between closely related species that are not manifested at the codon level.

Conclusions

Trpc4 was previously shown by a genetic screen to be a required gene for the regulation of NST by WSNs in mice. This study provided evidence from molecular evolutionary analyses that the regulatory regions of the TRPC4 protein evolved under episodic positive and purifying selection regimes in eutherian mammals consistent with selection on NST. These selection patterns implicate a discrete region of C-terminal sequence in the recruitment of TRPC4 into the regulation of NST, which should be a useful target of functional studies on NST. While TRPC4 may be a target of ongoing mammalian adaptation to diverse thermal niches, the data do not support a general increase in TRPC4 divergence rates among small-bodied mammalian clades that show greater phenotypic divergence in heterothermy than do larger-bodied outgroups.

Supplemental Information