Fourteen Carnivora species were chosen for evolutionary analysis based on the completeness of genome data, with humans as the outgroup (Figure 1). Initially, we obtained the full-length coding sequences (CDSs) of 39 HOX genes of 12 species, except for the red panda and northern fur seal, from the Orthologous Mammalian Markers database¹ (Douzery et al., 2014). Then, we examined the completeness of the start and stop codons, excluded potential pseudogenes, and eliminated sequences with any deletions or ambiguous sequences.

We obtained the hidden Markov model (HMM) corresponding to the HOX gene family (PF00046) (Finn et al., 2014) and used HMMER (version 3.1b2) (Prakash et al., 2017) for sequence alignment analysis of CDSs of the red panda. To identify homologous pairs of HOX genes, sequence alignment was performed using BLASTN (NCBI-blast-2.6.0+) (Chen et al., 2015) and BLASTP (NCBI-blast-2.6.0+) (Jacob et al., 2008). The database and query alignment included four rounds of mutual alignment. Walrus was used as the reference to identify the HOX genes of the northern fur seal. Similarly, to predict red panda HOX genes that were not annotated in the previous genome version, we performed homology prediction. We extracted the human HOX proteins from Ensembl (release 79) (Zerbino et al., 2018) and aligned them to the red panda genome using TBLASTN (version 2.2.23) (Gertz et al., 2006). Then, we extended the alignment regions by 10 kb at both ends and predicted the gene structure using GeneWise (version 2.2.0) (Birney et al., 2004). The optimal result was taken as the final HOX gene set of the red panda (Figure 2). Sequence analysis was carried out to check the integrity and lengths of all the predicted HOX genes.

Multiple alignments of the DNA and amino acid sequences of the HOX gene set were performed with two methods, MAFFT (version 7) (Katoh and Standley, 2013) and PRANK (Löytynoja, 2014) plus PAL2NAL (Suyama et al., 2006), to form the final gene set alignment. We defined the fourfold degenerate (4D) sites in the alignment sequence using custom scripts. The conserved sequence blocks were extracted from the outputs of multiple sequence alignment by Gblocks (version 0.91) (Castresana, 2000). The selection of the best substitution model for alignment was performed by ProtTest (version 3.4) (Darriba et al., 2011). The phylogenetic tree of the concatenated genes was generated in PhyML (Guindon et al., 2010). The amino acid substitution model was specified as GAMMA + JTTF. The obtained phylogenetic tree was visualized using FigTree viewer v1.4.1.²

To identify positively selected genes (PSGs), we conducted selective pressure analyses using CodeML from PAML 4.8 (Yang, 2007), with branch-site models employed (Yang and dos Reis, 2011). We divided our dataset into two large groups (G1 and G2): G1 included 10 species after removing four pinniped species and the sea otter and focused on the evolutionary analysis of HOX genes related to limb development of bamboo-eating giant and red pandas; G2 targeted the evolution of HOX genes associated with limb development in pinnipeds and sea otter and included 13 species after removing the giant and red pandas.

Furthermore, selective pressure analyses were performed using nine small groups under the same conditions: G1-a: setting the giant panda as the foreground branch; G1-b: the red panda; G1-c: both pandas; G2-a: the Weddell seal; G2-b: the Hawaiian monk seal; G2-c: the walrus; G2-d: the northern fur seal; G2-e: the Weddell seal, Hawaiian monk seal, walrus and northern fur seal; and G2-f: the sea otter. Two categories were used to set different foreground branches in pinnipeds (Figure 3): G2-e1, setting the most recent common ancestor of pinnipeds and the four pinniped species as the foreground branch; and G2-e2, setting only the four pinniped species as the foreground branch.

To confirm the positively selected genes identified in PAML analysis, we further used four models implemented in HyPhy from Datamonkey (Weaver et al., 2018) to examine the signatures of positive selection at three levels (site, gene and branch levels) based on the dN/dS ratio (ω). The analyses at the site level include three methods, which were used to identify individual sites subject to positive selection along subsets of phylogenetic tree branches: (1) FEL (Pond et al., 2005), (2) FUBAR (Murrell et al., 2013), and (3) MEME (Murrell et al., 2012). At the gene level, BUSTED analysis (Murrell et al., 2015) was used to identify genes for episodic diversification. At the branch level, the aBSREL method (Murrell et al., 2015) was used to identify positive selection on individual branches. Finally, to reduce false-positive results, only PSGs detected by both PAML and any one of the Datamonkey methods were considered.

We identified convergent amino acid substitutions between giant and red pandas and among marine mammals with the following rules: (i) the amino acid residues of both extant lineages were identical, and (ii) amino acid changes were inferred to have occurred between the extant lineage and its most recent ancestor lineage (Hu et al., 2017). In this study, calculation of frequencies and rates for the categories and reconstruction of ancestral protein sequences were performed by the CodeML program in PAML 4.8. The relative evolutionary rates of all sites within the gene followed the gamma distribution of sitewise rate variation and the frequency of all amino acids in each site of the gene. We also counted the convergent substitution events and calculated the convergence probability for each branch pair. To filter out noise resulting from chance amino acid substitutions, we performed a Poisson test to verify whether the observed number of convergent substitutions of each gene was significantly greater than the expected number caused by random substitution under the JTT-f gene and JTT-f site amino acid substitution models. We reconstructed the substitutions on all branches in a mammalian phylogeny containing the 15 mammals studied. We then tallied all convergent amino acid substitutions on the giant and red panda branches for all HOX genes, including all possible strictly convergent changes. To further check the strictness of the convergent amino acids, we took advantage of the 100 prealigned vertebrate genomes in the University of California Santa Cruz genome browser.³

Finally, we detected the common amino acid substitutions that considered only the terminal branches rather than the ancestral branches. Specifically, giant and red pandas have the same amino acids that differ from those of other species, and the Weddell seal, Hawaiian monk seal, walrus, northern fur seal and sea otter have the same amino acids that differ from those of other species.

To evaluate the signatures of rapid evolution on each HOX gene, we employed a branch model implemented in the CodeML module of PAML 4.8. Two different tests were conducted to identify lineage-specific effects in the evolution of clades: a two-ratio model versus a one-ratio model and a two-ratio model versus a free-ratio model. The overlap of significant genes from the two tests was thought to be under rapid evolution.

Negative selection plays an important role in maintaining the long-term stability of biological structure and function by purging deleterious mutations. We identified negatively selected genes according to the following criteria: dN/dS of the selected foreground branch is lower than dN/dS of all branches or the background branch, and the p value must be lower than 0.05.

To explore the structural effects of candidate sites, the predicted 3D structure of candidate HOX genes was generated using the trRosetta algorithm (Yang et al., 2020). To evaluate the impacts of amino acid changes on the overall protein structure, we calculated the predicted changes in folding free energy (ΔΔG) between ‘wild-type’ and ‘mutant-type’ amino acids using DynaMut (Yang et al., 2020). We predicted the effects of mutations on protein stability and flexibility. We obtained the PDB files from trRosetta and applied DynaMut on the online website.⁴ All of the protein 3D structure figures were generated using PyMOL.⁵

To characterize the functional impact of the mutant, we initially used Protein Variation Effect Analyzer (PROVEAN) (Choi and Chan, 2015) to predict the potential effect of an amino acid substitution. Second, we also obtained the 3D structure of particular proteins. The 3D structure of the wild-type protein was built in silico.

After a series of strict screenings and integrations, we obtained the HOX gene set from 14 Carnivora species with humans as the outgroup. All HOX gene sequences from the above species were intact, without premature stop codons or frame-shift mutations, indicating the presence of functional HOX proteins. The number of HOX genes varied among the 15 species (Supplementary Table S1), and not all 15 species included all members of the HOX gene family (Figure 1). These genes belong to four HOX gene clusters, each of which contains different numbers of HOX genes. Specific information on HOX genes, such as start codons, stop codons and gene lengths, is available in Supplementary Table S1.

For the red panda, only 25 genes could be downloaded from the previously published genome (Hu et al., 2017; Supplementary Table S1). We annotated an additional 14 HOX genes from the genome of the red panda. The results of sequence alignment showed high identities (average > 90%) to corresponding genes in humans (Figure 2). In total, we obtained 39 HOX genes of the red panda for subsequent evolutionary analysis.

Based on the concatenated CDSs of HOX genes of 15 species, we constructed an ML phylogenetic tree to decipher the evolutionary relationships of Carnivora species. We constructed a neighbor-joining phylogenetic tree using 4Dtv sites (Figure 1), which was consistent with the phylogenetic tree constructed based on nuclear and mitochondrial genes (Flynn et al., 2005) or on large-scale integrated analyses (Nyakatura and Bininda-Emonds, 2012; Upham et al., 2019).

By manually checking the multiple sequence alignment, we removed potential false-positive selection sites detected by PAML. When the giant panda, red panda, or both pandas were used as the foreground branch (G1-a; G1-b; G1-c), positive selection was detected for HOXA3 with one positively selected site (176:G:0.811, ω = 280.90, p < 0.001), HOXB4 (174:S:0.55, ω = 999, p < 0.001), HOXA3 (150:G:0.976*, ω = 97.02, p < 0.001) and HOXD4 (38,G:0.822, ω = 178.24, p < 0.001) (Figure 3A; Table 1).

For marine Carnivora species, when the sea otter was the foreground branch, only HOXA6 was detected to be under positive selection with one positively selected site (102:S:0.902, ω = 999, p < 0.001) (Figure 3A; Table 1). When the Weddell seal, Hawaiian monk seal, walrus or northern fur seal were used as the foreground branch separately, no significant PSGs were detected. However, we found several potential PSGs that were marginally significant, including HOXA6 (102:S:0.871; 211:D:0.813, p = 0.061) and HOXD12 (154:L:0.703, p = 0.067) when setting the Weddell seal as the foreground branch and HOXD3 (111:Q:0.833, p = 0.059) when the northern fur seal was used as the foreground branch (Figure 3A; Table 1).

Interestingly, the results of the two categories for the evolution of pinniped species were disparate (Figure 3A). Two HOX genes were identified under positive selection, with no common genes in the two categories. In category 1 (G2-e1), HOXB1 was identified as a PSG with one positively selected site (270:P:0.844, ω = 289.12, p < 0.05) during the origin and evolution of the pinnipeds. In category 2 (G2-e2), HOXB6 was identified as a PSG, with two positively selected sites (92:G:0.795; 140:S:0.998**, ω = 110.09, p < 0.001) (Figure 3A, Table 1).

We further analyzed the data using an additional method implemented in HyPhy of Datamonkey (Supplementary Figure S1 and Supplementary Table S2). When the giant panda was set as the foreground branch, three, seven and nine genes were detected at the gene, branch, and site levels, respectively. After integration, five genes, including eight sites, were detected under positive selection (Supplementary Tables S2–S4). In particular, HOXB3 was detected at all three levels. When the red panda was set as the foreground branch, we detected only one PSG (HOXB3), including two sites, at the branch and site levels (Supplementary Table S3). When both the giant and red pandas were set as the foreground branch, the integrated results showed that four PSGs were detected, including HOXA3, HOXB3, HOXD11, and HOXD12 (Supplementary Table S2). Within these genes, 7 sites were detected by at least two site-based methods (Table 1 and Supplementary Table S4). Codon 150 of HOXA3 was also detected in the PAML analysis.

When the sea otter was set as the foreground branch, HOXA6 was detected as a PSG by gene-level and branch-level HyPhy analyses (Supplementary Table S2) and in the PAML analysis (Table 1). When the Weddell seal was used as the foreground branch, two PSGs, HOXB6 and HOXC13, were detected by integrated methods (Supplementary Table S2). No PSGs were found when the Hawaiian monk seal, northern fur seal or walrus was set as the foreground branch.

When combining the pinniped species together, we detected common PSGs (HOXA6, HOXB6 and HOXC13) in either category 1 or category 2 (Supplementary Table S2). Additionally, codon 270 of HOXB1 was detected both in our site-based results and in PAML, and codon 140 of HOXB6 was found in the integrated analysis and in PAML (Table 1; Supplementary Tables S2–S4). All the above results suggested that a small number of HOX genes and sites were under positive selection in Carnivora.

To determine whether there was convergent evolution of HOX genes between the giant and red pandas and between the pinnipeds and sea otter, we identified the signatures of amino acid convergence in the above pairs. For the giant and red pandas, only one convergent amino acid change was detected for HOXC10. The observed number of convergence events (1) was significantly higher than the number of convergence events expected based on the phylogenetic distance (0.2771; p = 0.0307). The amino acid sequence alignment identified an amino acid substitution of Lys236Gln, with a codon change from AAA to CAA in the giant panda and from AAA to CAG in the red panda (Figure 4A). This substitution was also found to be unique to both pandas when not considering the amino acids of ancestral branches (Supplementary Table S5). In addition, in the gene alignment of 56 mammals, the amino acid at the 236 position is Gln only in the giant and red pandas, whereas this position is Lys in other species, highlighting the potential functional role (Figure 4A).

For the pinnipeds and sea otter, we did not detect convergent amino acid changes in any HOX gene. They did not have any common amino acids when only the extant lineages were compared. However, the four pinniped species had 13 common amino acids from 10 genes that distinguished them from other species, including the limb development-related genes HOXA10, HOXA13, HOXB9, HOXB13, HOXC10, and HOXD12 (Supplementary Table S5).

We detected rapid evolution and negative selection of HOX genes (Table 2; Supplementary Tables S6–S9). With the giant panda as the foreground branch, we identified two rapidly evolving genes, HOXA3 and HOXD4. With the red panda as the foreground branch, HOXD4 was shown to be under rapid evolution. When both giant and red pandas were used as the foreground branch, HOXA3 was detected to be rapidly evolving. When any one of four pinnipeds (Weddell seal, Hawaiian monk seal, walrus, northern fur seal) or sea otter was set as the foreground branch, no rapidly evolving genes were detected. When four pinniped species were used as the foreground branch simultaneously, HOXC13 was detected to be under rapid evolution (Supplementary Tables S6, S7).

We identified four genes under negative selection (Table 2; Supplementary Tables S8, S9), including HOXC12 for the giant panda, the red panda, both pandas, and the sea otter as the foreground branch; HOXA6 for both giant and red pandas; HOXD4 for four pinniped species, the walrus, and northern fur seal, respectively; and HOXA3 for the walrus.

For the positively selected gene HOXA3, the effect of the amino acid mutation was predicted to be neutral using either the giant panda or red panda protein sequence. The ‘neutral’ outcome was also predicted for HOXC10, a gene convergently evolved between giant and red pandas (Table 3). According to the Pfam database, the convergent amino acid substitution is not located in the Homeodomain. Furthermore, we constructed the 3D structure of HOXA3 and HOXC10 proteins using trRosetta modeling, a deep-learning method, and predicted the dynamic effects of the mutations. For HOXA3, the giant panda G338S mutation destabilized the HOXA3 protein by increasing the molecular flexibility (0.319 kcal/mol/k) and shifting the Gibbs free energy (ΔΔG) value to the negative range (−0.374 kcal/mol) (Figure 2B; Table 3). When using the red panda protein sequence, the G338A mutation decreased molecular flexibility (−0.547 kcal/mol/k) and stabilized the protein. For HOXC10, the K236Q mutation destabilized the protein (−0.085 kcal/mol) but decreased the molecular flexibility (−0.159 kcal/mol/k) (Figure 4B; Table 3).

HOXC10, a DNA-binding transcription activator, plays a key role in anterior/posterior pattern specification, proximal/distal pattern formation, embryonic limb morphogenesis, and skeletal system development. Similar to other HOX genes, the expression of HOXC10 has a very broad pattern in terms of temporal and spatial extent at different stages of embryonic development depending on the cell environment and internal state (Wellik and Capecchi, 2003). HOX9 and HOX10 function together to pattern forelimb stylopods, and HOX10 also affects certain phenotypes in zeugopods (Wellik and Capecchi, 2003). Based on Bgee (Bastian et al., 2021), HOXC10 is highly expressed in the triceps brachii, biceps brachii, cartilage tissue, trabecular bone tissue, and upper arm skin in humans.

The two pandas are not closely related, and their sharing of adaptive traits reflects convergent evolution. Sesamoid bones are small auxiliary bones that form near joints and contribute to their stability and function (Antón et al., 2006). Thus far, providing a comprehensive developmental model or classification system for this highly diverse group of bones has been challenging. Based on the latest research, sesamoid bones are regulated by both TGFβ and BMP signaling pathways, and both genetic and mechanical regulation are involved in facilitating developmental diversity (Eyal et al., 2019). All types of sesamoid bones originate from SOX9+/SCX+ progenitors under the regulation of TGFβ and are independent of mechanical stimuli from muscles (Eyal et al., 2019). Similarly, HOXC10 and SOX9 were found to be enriched in the knee at some overlapping time points (Pazin et al., 2012), and HOXC10 is involved in the TGFβ/BMP and Wnt signaling pathways. In addition, PITX1 also strongly associates with many functionally verified limb enhancers that exhibit similar levels of activity in the embryonic mesenchyme of forelimbs and hindlimbs (Park et al., 2014). PITX1 can induce the expression of TBX4, HOXC10 and HOXC11 in chick forelimbs and the expression of TBX4 and HOXC10 in mouse forelimbs (Logan and Tabin, 1999). TBX4, TBX5, and HOX cluster genes are crucial for forelimb development, and mutations in these genes are responsible for congenital limb defects (Jain et al., 2018). TBX4 and HOXC10 interact directly in limbs and synergistically activate transcription via a T-box–HOX composite DNA sequence, and the transcriptional activities of TBX4 and HOXC10 depend on their DNA-binding sites (Jain et al., 2018). This suggests that HOXC10 might play a role in the development of pseudothumbs. In addition, pinnipeds and manatees (belonging to the family Trichechidae) underwent parallel evolution of HOXC10 (Li et al., 2018), which implies a potential important function.

Furthermore, no direct studies have shown that HOXC10 interacts with PCNT or DYNC2H1, two genes that were identified as possibly related to pseudothumb development in a previous convergent evolution study of giant and red pandas (Hu et al., 2017). Thus, HOXC10 could be another candidate gene for pseudothumb development for future functional verification.

Generally, the rates of nonsynonymous and synonymous substitutions in HOX genes are relatively low due to the strong conservation of the HOX gene family. In this study, we identified only one gene, HOXA3, under positive selection in pandas. HOXA3 mutations may be related to parathyroid gland organogenesis and pharyngeal organ development (Gordon, 2018). Furthermore, rapid evolution analysis found that giant and red pandas had a common rapidly evolving gene, HOXD4. Previous studies suggested that HOXD4 acts in parallel to regulate the expression of target genes directing skeletogenesis (Folberg et al., 1999). HOXD4-transgenes are specifically activated in chondrocytes, and mutations in this gene cause severe cartilage defects due to delays in cartilage maturation (Kruger and Kappen, 2010).

For the negatively selected gene HOXC12, there have been few relevant studies. HOXC12 has undergone strong purifying selection, which suggests that mutations in this gene may be harmful to organisms. In addition, HOXA6 was detected to be under purifying selection when both pandas were considered as the foreground branch. Previous studies reported that HOXA6 is expressed at a high level in several types of malignant tumors (Dickson et al., 2013). Whether these genes are involved in body or limb development needs further study.

We found no signatures of positive selection or convergent amino acid substitutions of HOX genes between the sea otter and pinnipeds. This suggests that the phenotypic convergence of marine Carnivora species may be achieved through gene expression or regulatory region variations of HOX genes or the evolution of other relevant genes. Similarly, few signatures of common positive selection on HOX genes were detected across three marine mammalian lineages (pinnipeds, cetaceans, and sirenians), and convergence occurred at a functional level of HOX genes (Nery et al., 2016).

However, focusing only on the pinnipeds, HOXB1 was identified to be under positive selection during the origin and evolution of pinnipeds. This gene is part of a developmental regulatory system that provides cells with specific positional identities on the anterior–posterior axis (from the UniProt database). HOXB6 was identified as a positively selected gene when each pinniped species was set as the foreground branch. Regarding purifying selection, HOXC13 was detected as a negatively selected gene in the ancestral linage of four pinniped species. Studies have shown that HOXC13 has high expression in integument development, and its mutations are associated with skin and appendage development (Wu et al., 2013). Interestingly, this gene showed a rapid evolution signature in the order Sirenia (Wang et al., 2009). A cetacean study showed that the evolution of the cetacean forelimb may be associated with the positive selection or selective relaxation of HOXD12 and HOXD13 (Li et al., 2018).

Although selection pressure appears to have varied among different lineages, HOXD4 was under strong purifying selection in the northern fur seal and walrus lineages, even in the ancestral linage of the four pinniped species, suggesting that HOXD4 may have been under purifying selection over a long evolutionary time. HOXA3 was detected as a negatively selected gene for the northern fur seal. HOXC12 was subject to negative selection in the sea otter lineage and was also a negatively selected gene for both pandas, suggesting that its functional relaxation could be more strictly constrained.

In summary, our study explored the molecular evolution of HOX genes in Carnivora and focused on the potential relationship between HOX9 ~ 13 genes and limb development, providing insights into the potential molecular evolutionary mechanisms of Carnivora limb development. Overall, a few HOX genes undergo positive selection or convergent evolution, most likely because of the functional importance and evolutionary conservativeness of HOX genes. In our study, the identified PSGs and convergently evolved genes among HOX9 ~ 13 genes could be important candidate targets for further functional verification. A combination of evolutionary analyses and functional verification would illuminate the mechanisms of evolutionary developmental biology for specialized limbs or appendages (Hu et al., 2023).