For the analysis of the seven mammalian sirtuin protein family members SIRT1 to SIRT7 we obtained coding DNA sequences, including several isoforms or different prediction entries from NCBI Genbank (Clark et al., 2015). The sequences were retrieved for human (Homo sapiens), mouse (Mus musculus), rat (Rattus norvegicus), American beaver (Castor canadensis), and guinea pig (Cavia porcellus). The coding DNA sequence for SIRT5 for rats could not be retrieved. Additionally, we obtained some sirtuin coding DNA sequences for marmots from Genbank, including sequences of SIRT1, SIRT3, SIRT5, and SIRT6 for yellow-bellied marmot (Marmota flaviventris) and the sequence of SIRT 5 for Alpine marmot (Marmota marmota marmota, Gossmann et al., 2019). A complete list of retrieved sirtuin sequences along with their respective identifiers can be seen in Supplementary Table 1 (see Supplementary Material). To conduct an analysis tailored to marmots, additional publically available coding sequences for SIRT1 to SIRT7 homologs were predicted. For this, coding sequences were predicted from genomic data based on sequence homology to the available Refseq sirtuin marmot sequences using Exonerate v2.2.0 (Slater and Birney, 2005). The genomes of the woodchuck (Marmota monax), the Himalayan marmot (Marmota himalayana) and the Vancouver Island marmot (Marmota vancouverensis) were retrieved from Genbank. Furthermore, coding sequences for the Mongolian marmot (Marmota mongolia), long-tailed marmot (Marmota caudata) and Altai marmot (Marmota baibacina) could be obtained by using re-sequencing information (Bai et al., 2019) (http://www.marmotdb.org/). Specifically, coding information were obtained by using the vcf's information obtained from mapping short read data onto the Himalayan reference genome by Bai et al. (2019). Each species for which sequence data was used was included in the phylogenetic tree that was provided obtained from Timetree (Kumar et al., 2017).

To prepare the obtained sequences a customized python pipeline was designed based on available scripts from Afanasyeva et al. (2018). In brief, it contained general automated preparation steps, e.g., renaming data files, adjusting identifier information and storing the coding sequences for each sirtuin SIRT1 to SIRT7 from each species. Note, that this pipeline could be easily extended to handle large scale data, such as entire proteomes or across large trees (Afanasyeva et al., 2018; Yusuf et al., 2020). Then, each coding DNA sequence was translated into the corresponding protein sequence using the standard genetic code. All sequences containing a stop codon within the sequence were disregarded for further analysis. The translated sequences were then aligned using MUSCLE v3.8 using default parameter settings (Edgar, 2004). These alignments were used for a manual inspection and limited to the selection of one sequence per species and protein. This was necessary, because alignment artifacts may artificially create a signature of positive selection (Afanasyeva et al., 2018). Criterions for the exclusions of sequences included sequences with missing fragments in comparison to the human canonical sequence as well as sequences with potentially erroneous intron exon boundaries. The identifier of the selected sequences which were retrieved from Genbank are listed in Supplementary Table 2, while the selected sequences of the predicted sirtuin sequences for SIRT1 are shown in Supplementary Table 3 (Supplementary Material).

Alignments of the selected and translated sequences where further processed using PRANK v.170427 using default parameter settings (Loytynoja, 2013) and critical residues were masked with ZORRO (Wu et al., 2012) when the confidence score was smaller than nine. Combined with the previous coding DNA sequences the results were used as input for PAL2NAL v14 (Suyama et al., 2006) to create coding DNA sequence alignments. Thereby, output files with renamed identifier were created as input for PAML (Yang, 2007). Prepared codeml files and phylogenetic tree files pruned down to the number of available proteins were then used as input for codeml of the PAML package.

To conduct a selection analysis based on substitution rates, codeml from the PAML package version 4.9j was used. Two model types were considered. First, a branch model comparison was executed to determine heterogeneity in dN/dS between the marmot clade and the rest of the tree (two-ratio model) relative to a model where a single dN/dS ratio (one-ratio model) was assumed. Secondly, branch-site comparisons were used to test for the evidence of positive selection within the marmot clade. We compared the branch-site model A to a model where a potentially positively selected site category was restricted to drift (e.g., dN/dS = 1). Furthermore, the site models model 1 (M1a, NearlyNeutral, 2 categories) (Wong et al., 2004; Yang et al., 2005) and model 22 (M2a_rel, PositiveSelection, 3 categories) (Weadick and Chang, 2012), that allow the ω ratio to vary among sites (among codons or amino acids in the protein), were compared to a Clade model C (Bielawski and Yang, 2004).

The codeml-generated output files were further processed using a customized python script. For the branch test comparison dN/dS values were extracted from the one ratio and two ratio models and twice the log likelihood difference was assumed to be χ² distributed with one degree of freedom. Same accounts for the three different branch-site model comparisons. Here, the dN/dS values were extracted for model A, model A fixed, model=1 (M1a) and model=22 (M2a_rel), as well as for Clade model C. Again, twice the log likelihood difference of model A against model A fixed (with ω = 1) was assumed to be χ² distributed with one degree of freedom. For model M1a against Clade model C d.f.=3 was used while one degree of freedom was used for model=22 (M2a_rel) against Clade model C. Furthermore, the NEB (Naïve Empirical Bayes) and BEB (Bayes Empirical Bayes) values were extracted from rst-files from the codeml output for data plotting.

A 3D model was calculated using the SWISS-MODEL (Guex et al., 2009) (with access on SWISS-MODEL Interactive Workspace, n.d.) a template for human SIRT1 protein structure from Protein Data Bank (4I5I) (Zhao et al., 2013) and the coding sequence of human SIRT1. With FirstGlance (FirstGlance in Jmol, n.d.) the model was visually processed. The influence of the substitution of selected amino acids regarding the secondary structure of the protein was analyzed using Jpred4 (Drozdetskiy et al., 2015) as well as Porter 5.0 (Torrisi et al., 2019) submitting the protein sequence of human SIRT1 in comparison to the same sequence with single amino acid substitution [with substitution of lysin (K) to threonine (T) in amino acid sites 499, serine (S) to leucin (L) in amino acid sites 508 and substitution of arginine (R) to isoleucine (I) in amino acid sites 613] deducted from the marmot sequences of M. sibirica, M. himalayana and M. caudata, and M. flaviventris, respectively.

Differences in the evolutionary rates, recovered from between-species comparisons, may reveal molecular footprints of selective pressures. Here we investigated heterogeneity in the rate of molecular evolution of sirtuins between the marmots and other mammalian species (Figure 1).

The ratio of non-synonymous to synonymous substitutions ω = dN/dS was used as a measure for the rate of molecular evolution at the protein level. The dN/dS ratio indicates the driving selective forces in the evolutionary process, as elevated dN/dS ratios above one suggest the influence of positive selection (ω > 1) while a low ratio (ω <1) reveals a process driven mainly by negative selection (Gossmann and Schmid, 2011). The comparison of the estimated evolutionary rate between the one ratio (single dN/dS ratio for all considered species) and two ratio model (different dN/dS ratio for marmot clade) revealed higher dN/dS ratios in the marmot clade relative to the compared species for six of the seven sirtuin proteins (SIRT1, SIRT3, SIRT4, SIRT5, SIRT6, and SIRT7) (Table 1).

Indeed, we find a significant difference of the d_N/d_S ratios between the marmot clade and the remaining tree for SIRT1. This shows that the marmot clade has a different evolutionary rate regarding SIRT1 than the remaining branches of the phylogenetic tree. The explanation for the higher dN/dS ratios could be relaxed purifying selection, but it could also suggest a possible influence of positive selection specific to marmots.

Since the branch test revealed an increased rate of molecular evolution in the marmot clade relative to other mammalian species, a possible role of positive selection events was further analyzed. To investigate positive selection in the evolution of SIRT1 in the marmot clade three branch-site model comparisons were executed. The branch-site model A was compared to the model A fixed where a potentially positively selected site category was restricted to drift (with dN/dS = 1). Furthermore, Clade C was compared to model 1 (M1a, NearlyNeutral) and also to model 22 (M2a_rel, PositiveSelection). All three comparisons revealed a significant difference (Figure 2), indicating the occurrence of positive selection events within the marmot clade.

A significant difference for the three different likelihood comparisons is assumed if twice the log likelihood difference (2Δl = 2 · (l1 – l2) with l1 and l2 as the log of the maximum likelihood (ML) estimated of two models) is χ² distributed. A likelihood comparisons between the branch-site model A and model A fixed, where a potentially positively selected site category was restricted to drift (e.g., dN/dS = 1) was performed to test for the evidence of positive selection within the marmot clade (χ² distribution with d.f. = 1). Furthermore, results of a likelihood comparisons between Clade model C (Bielawski and Yang, 2004) and site models model 1 (M1a, NearlyNeutral, 2 categories) (Wong et al., 2004; Yang et al., 2005) (χ² distribution with d.f. = 3) and a comparison between Clade model C and model 22 (M2a_rel, PositiveSelection (3 categories)) (Weadick and Chang, 2012) (χ² distribution with d.f. = 1) differed significantly.

To identify individual amino acid substitutions that were likely promoted by positive selection, we determined every amino acid's probability to belong to a class with either low (dN/dS <1), neutral (dN/dS = 1) or a high dN/dS ratio (dN/dS > 1). As the different models revealed similar results, we focused on model A results. The gained NEB (Naive Emperical Bayes) probabilities show three amino acid sites with high posterior probabilities (> 75%) for site class 2 with ω = dN/dS = 1.68 (Figure 3), indicating that the footprints of positive selection can be pinpointed to particular amino-acid replacements. Similar results have been obtained for posterior probabilities based on BEB (Bayes empirical Bayes, Figure 3) probabilities. This approach accounts for sampling error (Yang et al., 2005) illustrating that the amino acids sites underlying the signal of positive selection can be robustly identified.

The three outlined amino acid sites presumably have undergone positive selection. The alignment of the analyzed sequences for SIRT1 shows the amino acid sites and the corresponding substitutions within the marmot clade (Figure 4). In respective to coding sequence of human SIRT1 the amino acid sites are 499 with substitution of lysin (K) to threonine (T) in M. sibirica, 508 with substitution of serine (S) to leucin (L) in M. himalayana and M. caudata and 613 with substitution of arginine (R) to isoleucine (I) in M. flaviventris.

Since we identified substitutions that presumably evolved under positive selection, we were interested whether we could identify a possible impact of these residues on the structural level. To investigate this effect further, a secondary structure prediction of human SIRT1 protein was performed along with in-silico mutated variants of the human protein, substituting the three residues we identified from the marmot species. Intriguingly, we find that a differently predicted substructure within the amino acid sites 518 to 523 of the human SIRT1 sequence is observed for each of the 3 residues (Figure 5). Interestingly the secondary structure of this 6 amino acid long sequence stretch can be predicted only with low confidence, but it is embedded in a coiled coiled region that can be confidently predicted. The observed effect on the secondary structure is also remarkable, because one of the mutations is 270 bp (90 aa) away from the affected amino acid stretch. While the human sequence is predicted with a helix structure for the mentioned amino acid sites, the sequences with the single substitutions each are predicted with a predominantly coiled structure instead. We conducted control mutations where we randomly substituted amino acid residues nearby the identified positions and did not find such an effect on the structure prediction for serine and arginine (Supplementary Figures 3E–G). We find that changes of nearby lysin residues to the observed mutations had a similar impact that lead to a coiled-coiled prediction in the target region (Supplementary Figures 3H–J).

This structural difference may alter the enzyme function. As Zhao et al. (2013) reported the catalytic domain of SIRT1 reaches up to the residue 516. Their estimated structure (shown in Supplementary Figure 1) which ends with the amino acid E 512 suggests a very close position of the subsequence QKELAY with amino acid sites 518 to 523 to the substrate bounded in the catalytic domain. Assuming a coiled structure for the subsequence (as predicted in Figure 5) instead of the shown strand in Supplementary Figure 1 it is conceivable that the subsequence QKELAY might have an influence on the substrate affinity.