Transcript sequences encoding the longest peptide for established MMPs from representative tetrapod taxa (human, mouse, Xenopus, chicken) and a tetrapod outgroup (Latimeria) were retrieved from NCBI¹ and (Fu et al., 2009), who previously characterized MMPs from X. tropicalis. We used these sequences as queries in Blast searches (tBLASTX) (Altschul et al., 1990) against publicly available axolotl genomic and transcriptomic databases and datasets. These included Sal-Site² (Smith et al., 2005; Baddar et al., 2015) and axolotl-omics³ (Nowoshilow and Tanaka, 2020), in addition to transcriptome datasets from Bryant et al. (2017) (GSE92429), (Dwaraka et al., 2018) (GSE116615 and GSE116777), (Nowoshilow et al., 2018) (PRJNA378970, PRJNA378982), and (Smith et al., 2019) (GCA_002915635.2), and (Schloissnig et al., 2021) (PRJNA520877, PRJNA644663, and PRJNA645452). Axolotl sequences that yielded significant alignments (query coverage > 90%, highest bit score) were subsequently used to run additional blast searches against the aforementioned axolotl transcriptomic/genomic databases using BLASTN (query coverage > 90%) to identify potential duplicates within axolotl. Presumptive axolotl MMP transcript sequences were then manually curated to remove splice variants and transcripts encoding partial protein sequences. This list of full-length axolotl MMP transcripts was then used to search for homologs in other salamander species. This was performed by using the NCBI tBLASTN search program (query coverage > 90%, highest bit score) against available transcriptome assemblies of the following salamanders: [Ambystoma texanum, Ambystoma laterale, Ambystoma tigriunum (McElroy et al., 2017)], and Hynobius chinensis (Che et al., 2014), and by using Blast2go software (Götz et al., 2008) to search BLASTN/BLASTX transcriptome assemblies of other salamanders: Ambystoma maculatum (Burns et al., 2017; Dwaraka et al., 2018), Ambystoma andersoni (Dwaraka et al., 2018), Cynops pyrrhogaster (IMORI)⁴ described in Casco-Robles et al. (2018), Pleurodeles waltl (Elewa et al., 2017), Nothophthalmus viridescens (Abdullayev et al., 2013; Elewa et al., 2017), and Bolitoglossa ramosi (Arenas Gómez et al., 2018). Full-length transcripts were included in the analyses and partial sequences were discarded. To gain additional amphibian phylogenetic perspective, we used axolotl MMP presumptive sequences as queries to identify MMP orthologs in three caecilians: (Rhinatrema bivittatum, Typhlonectes compressicaud, and Microcaecilia unicolor) (Torres-Sánchez et al., 2019; and NCBI, see text footnote 1). Overall, a total of 268 MMP full coding sequences were used and translated into amino acid sequences using the NCBI ORF tool.⁵ ClustalW (Madeira et al., 2019) was used to align sequences and annotate protein domains using human MMPs as protein models. Source information identifiers for the MMPs discovered in this study, including established genome locations for axolotl MMPs (Nowoshilow et al., 2018; Smith et al., 2019), are shown in Supplementary Tables 1, 2. This same general approach was also used to extract 60 TIMP sequences from the axolotl and all the taxa above (Supplementary Table 3). Axolotl MMP and TIMP gene annotations were determined following the guidelines recently described in Nowoshilow et al. (2021).

Protein sequences of the 268 MMPs were aligned using ClustalW (Madeira et al., 2019) to generate a multiple sequence alignment (MSA) using RevTrans (Wernersson and Pederson, 2003). The quality of the resulting MSA was examined using AliView to ensure alignment of conserved motifs across all the MMP sequences (Larsson, 2014). Phylogenetic analyses were conducted using the IQ-tree command line tool to construct a maximum likelihood tree under a GTR+F+R8 substitution model identified by IQtree built-in ModelFinder with 100,000 bootstrap support value replicates (Nguyen et al., 2015; Kalyaanamoorthy et al., 2017; Hoang et al., 2018). The phylogenetic tree was visualized and annotated using Figtree V1.4.2 software⁶. This same approach was used to construct a TIMP sequence evolutionary tree under a TN+F+I+G4 substitution model identified by IQtree built-in ModelFinder with 100,000 bootstrap support value replicates (Nguyen et al., 2015; Kalyaanamoorthy et al., 2017; Hoang et al., 2018). 3D structural alignments of human and axolotl sequences were generated using PyMOL (The PyMOL Molecular Graphics System, Version 1.2r3pre, Schrödinger, LLC).

Axolotl BAC clones were isolated from existing libraries (Smith et al., 2009; Voss et al., 2013) and sequenced to verify the identity and genomic locations of 5 novel axolotl MMPs. Briefly, axolotl BAC clone superpools were screened using PCR primers (Supplementary Table 4) that were designed from Sal-Site EST contigs (Smith et al., 2005; Baddar et al., 2015) to generate amplicons for three MMP3/10 and two MMP13 paralogs. The location of positive clones among BAC library microtiter plates was determined by sequential PCR of plate, column, and row BAC pools. BAC clone DNA was isolated using the PureLink HI Pure Plasmid Maxiprep Kit (Invitrogen) and then sequenced on a PacBio RSII by the Duke Center for Genomic and Computational Biology. Sequences were analyzed using DNASTAR SeqMan (DNASTAR, Inc., Madison, United States).

For MMP gene expression analysis, RNA seq reads were mapped to the recent release of the axolotl genome AmexG_v6.0-DD (Schloissnig et al., 2021) using HISAT2 aligner v.2.2.0 (Kim et al., 2019). Depth of coverage normalized by reads per million was computed with bedtools v2.27.1 (Quinlan and Hall, 2010). Average FKPM values were calculated for each located MMP across three replicas of datasets corresponding to wound epidermis (SRR7499357, SRR7499358, SRR7499359), two replicas of distal blastema (SRR2885553, SRR2885591) and two replicas of proximal blastema (SRR2885865, SRR2885866).

A comprehensive survey of MMP genes from available axolotl transcriptomic, genomic and EST databases, along with targeted sequencing of presumptive axolotl novel MMPs genes, yielded a total of 28 MMPs (Supplementary Table 1). This number exceeds the number of MMPs identified for other tetrapods, including anciently duplicated teleost genomes (26 MMPs are known for zebrafish) (Pedersen et al., 2015). This survey shows that the axolotl has an expanded number of MMP paralogs, the highest number recorded for any species.

We next examined the evolutionary history of MMPs using 268 complete MMP protein coding sequences from 19 vertebrate species: (A) 10 salamander species: the axolotl (Ambystoma mexicanum (Amex), Ambystoma andersoni (Aand), Ambystoma texanum (Atex), Ambystoma laterale (Alat), Ambystoma tigrinum (Atig), Hynobius chinensis (Hchi), Pleurodeles waltl (Pwal), Nothophthalmus viridescens (Nvir), Bolitoglossa ramose (Bram), Cynops pyrrhogaster (Cpyr), (B) three caecilians [Rhinatrema bivittatum (Rbiv), Typhlonectes compressicaud (Tcom), Microcaecilia unicolor (Muni)], (C) Xenopus tropicalis (Xtro), (D) Gallus gallus (Ggal), (E) Mus musculus (Mmus), (F) Homo sapiens (Hsap), and (G) Latimeria chalumnae (Lcha) as a tetrapod outgroup (Figure 2 and Supplementary Table 1). The majority of clades show overall strong bootstrap values indicating that the inferred relationships within and between the MMP subfamilies are well-supported. More than one MMP sequence was identified for salamander species within different clades, consistent with salamander specific gene duplication. These clades include: salamander-specific collagenases, salamander-specific stromelysins, MMP13, and novel MMPe. We discuss each of these clades below.

Archetypal MMPs share the same protein domain architecture (Figure 1) and can be divided into three subcategories: collagenases (MMP1, 8, 13), stromelysins (MMP3, 10), and others (MMP12, 19, 20, 27) (Nagase et al., 2006; Page-McCaw et al., 2007; Jackson et al., 2010). The distribution of archetypal MMPs was heterogeneous in the MMP tree (Figure 2). For example, MMP13, MMP19, and MMP20 clades are separate from each other and each is inclusive of all corresponding Latermeria/tetrapod orthologs. The remaining archetypal MMPs grouped into several amphibian and non-amphibian clades. Specifically, the clade annotated as “Non-amphibian collagenases and stromelysins” included the following collagenases and stromelysins from non-amphibian tetrapods: MMP1 orthologs in human and chicken, murine specific MMP1 duplicates (MMP1A and MMP1B), mammalian MMP8, MMP12 orthologs, and MMP3, MMP10, and MMP27 found in mouse, human, and chicken. A separate clade annotated as “amphibian novel collagenases” (Figure 2) included some amphibian collagenases from Xenopus, caecilians, and salamanders. This clade comprises a small subset of MMP1 orthologs in Xenopus and caecilians, a novel collagenase isolated in Xenopus (MMP18) (Stolow et al., 1996) and its likely ortholog in caecilians (MMP18l), and a salamander-specific collagenase subclade with four different collagenases: COL (A–C). A previously cloned newt-collagenase (nCOL, GenBank: AAX14806; Vinarsky et al., 2005; Supplementary Table 1) is predicted to be an ortholog of Amex_COL A. The identification of COL A–C orthologs for more than one salamander species suggests these are distinct loci and we note that loci encoding COL A, B, and C are syntenic on axolotl chromosome 7, consistent with tandem duplication (Supplementary Table 2).

The clade of the third collagenase (MMP13) also shows evidence of gene duplication within amphibian lineages. While Xenopus and other vertebrate species have a single MMP13, caecilians and salamanders seem to encode an additional paralog. The clade topology suggests that MMP13 paralogs are ancestral to caecilians and salamanders. We identified and sequenced the genomic sequences of Amex_MMP13A and Amex_MMP13B. Both genes are tandemly located on the same chromosome (Figure 3 and Supplementary Table 2). All salamander collagenases have domains characteristic of other tetrapod collagenases (Supplementary File 1), thus validating these gene name annotations.

Humans encode three stromelysins: MMP3, MMP10, and MMP11. Of these three stromelysins, MMP11 is distinguished by having a furin-recognition domain (Figure 1). While only a single MMP11 was identified among salamander species (Figure 2 and Supplementary Table 1), many gene duplicates were identified for the other two stromelysins. Until our study, few stromelysins (MMP3/10) were known for newts and salamanders (Miyazaki et al., 1996; Vinarsky et al., 2005; Supplementary Table 1 and Figure 2). Our findings suggest salamanders uniquely encode five stromelysins (MMP3/10 A–E) that share all characteristic protein domains in typical tetrapod stromelysins (MMP3 and 10) (Supplementary File 2). However, they are clustered in a separate clade from non-amphibian stromelysins, highlighting significant sequence divergence. Two previously identified stromelysins from C. pyrrhogaster- MMP3/10a and MMP3/10b (Miyazaki et al., 1996) (GenBank: D82053.1, D82054.1) and from newt N. viridescens (GenBank: AAX14804.1, AY857754.1) (Vinarsky et al., 2005) are likely orthologs of axolotl MMP3/10 A and B (Supplementary Table 5). Additionally, we isolated and obtained the genomic sequences of MMP3/10 B, C, and D from axolotl BAC clones to show that these are distinct loci in the axolotl genome. We could not identify orthologs for the remaining axolotl stromelysins (MMP3/10 C, D) in other salamanders, perhaps reflecting incomplete transcriptome sequencing. Four out of five axolotl MMP3/10 duplicates are tandemly located on chromosome 7 indicating they likely arose via tandem duplication (Supplementary Table 2).

In 2003, a novel MMP (nMMPe) was identified from C. pyrrhogaster regenerating limbs (Kato et al., 2003). This novel gene encodes a 502 amino acid protein and lacks homology with other vertebrate MMPs. Our search identified candidate nMMPe orthologs in the axolotl, H. chinensis, and A. andersoni (Figure 2 and Supplementary Table 1). We note that all of these novel MMPs cluster together and form a separate clade. MMPe lacks transmembrane insertions, suggesting it may belong to the archetypal secreted MMP category and not the cell-membrane anchored. This novel MMP most likely arose uniquely in the salamander lineage.

We note that while human and other non-amphibian tetrapods encode a single gene for each MMP17 and MMP21, salamanders, Xenopus, and caecilians encode an extra copy of each (Figure 2, clades: MMP17, MMP21). This suggests that perhaps a-specific duplication of these genes occurred at the base of amphibian ancestors. Alternatively, these gene duplicates were common to all tetrapods but were lost after amphibians’ divergence.

As shown in Figure 2, salamanders have orthologs to established MMPs in other tetrapods and they encode the protein domains characteristic of each: [MMP2, MMP9 (Supplementary File 3), MMP7 (Supplementary File 4), MMP11, 21l, 28 (Supplementary File 5), MMP14, 15, 16, 24 (Supplementary File 6), MMP17, 25 (Supplementary File 7), and MMP23B (Supplementary File 8)]. Most of these clades contain a presumptive ortholog from one or more of the representative species used in this study. No salamander-specific expansion is seen in any of these aforementioned clades. This pattern suggests a more conservative evolutionary history for these MMPs. Salamanders seem to lack orthologs for MMP8, MMP12, and MMP27. Because Xenopus and caecilians also lack orthologs for these genes, it seems likely that they were not present in the common ancestor of amphibians.

In summary axolotls and other salamanders are predicted to encode 28 MMPs (Supplementary Table 5), 9 of which are novel paralogs (collagenases, stromelysins) that arose uniquely in the salamander lineage. More importantly, all the novel salamander MMPs are of the secreted type, rather than cell membrane anchored. This indicates that salamanders evolved an extensive battery of secreted type MMPs that presumably execute catalytic activities in the extracellular milieu.

We found that salamander novel MMPs (MMP3/10 A-E, MMPe) have unique salamander-specific insertions and substitutions not found in other tetrapod MMPs. Unlike collagenases, stromelysins have longer hinge domains that play a role in substrate specificity that distinguishes these two subclasses of MMPs (Fu et al., 2009; Manka et al., 2019). The hinge domain length in all five axolotl collagenases (COL A–C, MMP13 A-B) is highly similar to that in other tetrapod collagenases, whereas all presumptive salamander stromelysins (MMP3/10 A–E) have variable insertions confirming their assignment as stromelysins (Supplementary File 9). Given the hinge domain role in substrate specificity in stromelysins, it is possible that salamander novel stromelysins bind different classes of substrates.

Salamander novel MMPe proteins also have unconventional variation in the conserved cysteine switch and hinge domain (Supplementary File 10). The conserved cysteine switch domain (PRCGVPD) is important for enzyme latency and residues surrounding it stabilize the interaction and therefore are highly conserved in most MMPs (Sanchez-Lopez et al., 1988; Van Wart and Birkedal-Hansen, 1990). Mutations in some of these surrounding residues were shown to weaken enzyme latency. For example, the substitution of the first proline by leucine results in activation of the proenzyme (Sanchez-Lopez et al., 1988). Salamander MMPe genes share a serine residue replacing the first proline in the cysteine switch (SRCGVPD). The effect of this substitution on MMPe function remains to be elucidated. The second distinctive feature of MMPe is the length of the hinge domain (Supplementary File 10). MMPe enzymes have the longest threonine- rich insertions in the hinge domain among all archetypal tetrapod MMPs. As hinge domains are usually critical for substrate binding and MMP catalytic efficiency, these novel unconventional sequence signatures may affect enzymatic activity in ways that are unique to salamanders.

Orthologous MMPs in salamanders also seem to have evolved novel insertions. We found that salamander MMP16 encodes a unique 27 amino acid long insertion in the catalytic domain close to the canonical Zn⁺² binding site (Supplementary Files 6, 11). This insertion is salamander-specific and not found in other tetrapods. MMP16 is membrane-anchored (Figure 1) and this 27 amino acid-insertion could potentially affect catalysis and substrate recognition of the enzyme.

Human MMP1 binds collagen II at 6 different sites; three at the catalytic domain, one at the hinge domain, and two in the hemopexin domain (Bertini et al., 2012). This binding initiates conformational changes, facilitated by the hinge domain, to allow subsequent hydrolysis of collagen by the activity of the active site (Fasciglione et al., 2012). The axolotl 4 collagenases (COL A–C) shared highest sequence homology (∼56%) with human MMP1 (Figure 4), so we were curious to find how conserved these collagen binding domains are among salamander novel collagenases. The predicted 3D structure of the aligned human MMP1 and the four axolotl collagenases (a–c) is shown in Figure 5, where purple and green shaded amino acids represent conserved and variable residues, respectively. The overall pattern of alignment indicates amino acid variation among multiple sites between human MMP1 and the axolotl collagenases. Variation at site 3, which is located closely to the Zn⁺² binding domain, suggests potential differences in the nature of unwinding and hydrolyzing collagen. This variation may mirror evolutionary differences in salamander collagen amino acid sequences and/or 3D conformation in both systems, as the specific sites on the enzyme have to fit specific amino acid sequences in collagen to dock the enzyme onto the substrates for subsequent hydrolysis.

We sought to compare the level of gene expression patterns between expanded and non-expanded salamander MMPs. To accomplish this, we extracted gene expression profiles for all available mmps during limb regeneration using normalized datasets corresponding to wound epidermis and blastema (Bryant et al., 2017; Dwaraka et al., 2018; Supplementary Tables 6, 7), and single cell RNA seq (sc RNA-seq) data (Gerber et al., 2018; Leigh et al., 2018; Rodgers et al., 2020; Li et al., 2021). As summarized in Supplementary Tables 6, 7, most of salamander expanded MMPs: [COl (A-B), MMP13 (A, B), MMPe, MMP3/10 (A–D)] are highly upregulated after limb amputation in the wound. Some are expressed in immune cells, including macrophages and neutrophils that are resident or recruited to the wound site: COL (A-B) (Miyazaki et al., 1996; Kato et al., 2003; Vinarsky et al., 2005; Gerber et al., 2018; Leigh et al., 2018; Rodgers et al., 2020; Li et al., 2021). On the other hand, most of the salamander MMPs that belong to the non-expanded subfamilies (MMP 7, 11, 19, 14, 15, 16, 23, 24, 25, 28) exhibit relatively modest upregulation during wound healing, blastema formation, and blastema outgrowth with the exception of MMP9, 14. However, MMP7, 9, 19, and 17l are expressed in macrophages at 1 and 6 days post-amputation (Rodgers et al., 2020). Taken together, these results demonstrate that almost all salamander expanded MMPs are more actively engaged in cellular activities during critical early phases of limb regeneration than non-expanded MMPs. The role of the novel identified insertions in these salamander-specific MMPs may be important for their activities during regeneration.

MMPs can degrade a plethora of structural and non-structural ECM components (Sternlicht and Werb, 2001; Nagase et al., 2006; Huxley-Jones et al., 2007; Caley et al., 2015). MMPs may have duplicated to accommodate the increasing variety of ECM components that arose during the evolution of vertebrate tissue complexity (Fanjul-Fernández et al., 2010). According to this hypothesis, the novel expansion of MMPs in salamanders should correlate with an expansion of MMP substrates. Indeed, we found that multiple genes encoding matricellular proteins, known to be targets to MMP activities, seem to have duplicated in salamanders. Specifically, we identified gene expansions in two known families encoding matricellular proteins: CCN (cellular communication network factors) and SPARC (secreted protein acidic and rich in cysteine), and an extra coagulation factor (F10b) (Supplementary Table 8).

In mammals, the CCN gene family encodes six secreted matricellular proteins involved in several biological processes including wound healing, cell migration, mitogenesis, adhesion, and ECM remodeling (Lipson et al., 2012; Krupska et al., 2015). These include: cysteine rich 61 (CYR61), connective tissue growth factor (CTGF), nephroblastoma overexpressed (NOV), and three Wnt-inducible secreted proteins (WISP1, WISP2, and WISP3) (Krupska et al., 2015). Structurally, all encoded CCN members are comprised of four cysteine-rich protein modules: (N-terminal signaling peptide, an insulin-like growth factor binding protein (IGFBP), a Willebrand type C repeat (VWC), a thrombospondin type 1 domain (TSP-1), and a cysteine knot carboxyl terminal (CT)) (Holbourn et al., 2008; Krupska et al., 2015). These domains are separated by linker regions that are prone to proteolytic cleavage by MMPs and other enzymes. Consequently, the liberated domains act as protein modules with pleiotropic biological roles (Holbourn et al., 2009). Salamanders have novel additional members of this gene family. The first was discovered in regenerating hearts of N. viridescens (newt-specific CCN, nsCCN) (Looso et al., 2012). Although we could not identify a true ortholog for this ns-CCN in the axolotl, we found another transcript (nsCCN-homolog) that shares relatively lower protein sequence similarity (<65%). This axolotl sequence has a novel insertion in the vWC domain not found in other mammalian CCN members (Supplementary File 12). In addition, we found that salamanders encode a second CTGF (CTGF-b) and CYR61 (CYR61-like) gene (Crowner et al., 2019), raising the total number of CCN gene family members to 10.

Similar to CCN protein, SPARC (secreted protein acidic and rich in cysteine) family members are composed of modules separated by linker regions that are cleaved by MMPs, which leads to their activation. Human MMP3, for instance, can cleave and thereby activate SPARC domains (Manka et al., 2019). A member of this gene family, extracellular matrix protein 2 (ECM2), was first identified in human (Nishiu et al., 1998) and encodes a 699 amino acids long peptide comprised of the following domains: [signaling peptide, integrin-binding sequence, a von Willebrand factor (vWFC), and the leucine-rich-repeat domain]. We identified another ECM2 transcript (ECM2-b) in the axolotl -and other salamanders. Axolotl ECM2 and ECM2-b encode 714 and 720 amino acid long peptides, respectively, but share relatively low protein sequence identity (41%) indicative of significant sequence divergence. Notably, other amphibians are also predicted to encode an ortholog for ECM2-b. For example, the Nanorana parkeri ECM2-like (XP_018425624.1) protein shares higher sequence identity with axolotl ECM2-b. This suggests that Ecm2 may have duplicated prior to the divergence of salamanders and anurans.

Almost all MMPs can dissolve fibrin clots and therefore facilitate wound healing (Sternlicht and Werb, 2001; Monaco et al., 2007). Injuries (extrinsic and intrinsic) stimulate formation of a fibrin clot by activating a cascade of clotting factors. The human genome encodes 12 coagulation factors (F1–F12) and F10 represents an important shared point where both extrinsic and intrinsic coagulation pathways converge (Smith et al., 2015). Once activated, F10a can ultimately activate prothrombin (F2) to thrombin (F2a) which in turn converts soluble fibrinogen to an insoluble fibrin clot. We found two F10 genes (F10 and F10-b) in the axolotl and other salamanders (Supplementary Table 8). The functions of these genes and novel MMPs may intersect to regulate clot formation and dissolution.

Tissue Inhibitors of MMPs (TIMPs) are known to inhibit the activities of MMPs (Nagase et al., 2006; Huxley-Jones et al., 2007). Four different TIMPs are known for vertebrates (1–4). Given expansion of salamander MMPs, we explored the possibility that TIMPs might have expanded in parallel. On the contrary, only 5 TIMPs (TIMP1, TIMP1-like, TIMP2, TIMP3, TIMP4) are found in the axolotl genome (Figure 6), all of which have highly homologous orthologs in other tetrapods. Thus, the expansion of the salamander MMP repertoire is not coincident with a co-evolutionary expansion of TIMPs.