The axolotl genome and the evolution of key tissue formation regulators

Salamanders serve as important tetrapod models for developmental, regeneration and evolutionary studies. An extensive molecular toolkit makes the Mexican axolotl (Ambystoma mexicanum) a key representative salamander for molecular investigations. Here we report the sequencing and assembly of the 32-gigabase-pair axolotl genome using an approach that combined long-read sequencing, optical mapping and development of a new genome assembler (MARVEL). We observed a size expansion of introns and intergenic regions, largely attributable to multiplication of long terminal repeat retroelements. We provide evidence that intron size in developmental genes is under constraint and that species-restricted genes may contribute to limb regeneration. The axolotl genome assembly does not contain the essential developmental gene Pax3. However, mutation of the axolotl Pax3 paralogue Pax7 resulted in an axolotl phenotype that was similar to those seen in Pax3−/− and Pax7−/− mutant mice. The axolotl genome provides a rich biological resource for developmental and evolutionary studies.


A long-read assembler for large genomes
Our aim was to generate a genome sequence assembly for the d/d axolotl strain (Fig. 1a), which is commonly used in laboratory regene ration studies owing to its compatibility with live imaging. Taking into consideration the expected challenge of assembling the complex 32-Gb genome 6 , we sequenced 110 million long reads (32× coverage, N50 read length 14.2 kb) using Pacific Biosciences (PacBio) instruments (Supplementary Information section 1) to avoid the read sampling bias that is often found when using other technologies and to span repeat-rich genomic regions that cause breaks in short-read assemblies (Fig. 1b, c).
We developed an assembly algorithm (MARVEL) that integrates a two-phase read-correction procedure that keeps long PacBio reads intact for assembly ( Supplementary Information section 2). MARVEL produced a contig assembly with an N50 of 218 kb. Next, we used 7× Illumina-based sequencing to correct sequence errors in 1% of the contig bases (Fig. 1b), which yielded a sequence accuracy of more than 99.2%. On the basis of the Illumina data, we estimated a heterozygosity of 0.47% (Supplementary Information section 2.2).
To provide a scaffold for the contig assembly, we generated de novo optical maps using the Bionano Saphyr system (Supplementary Information section 2.3). The Bionano map resolved contig chimaeras, which were found in 1.7% of contigs, slightly reducing N50 contig length to 216 kb (Fig. 1d). The final hybrid assembly yielded an N50 scaffold length of 3 Mb. Compared to the short-read assembly of the 20-Gb spruce genome 7 or the 22-Gb loblolly pine genome 8 , which involved 12× long-read coverage, the axolotl assembly showed 56-and 29-fold improvements in contiguity, respectively (Table 1).
To assess the completeness of the assembly (Supplementary Information section 4.1), we first determined the number of aligning non-exonic ultraconserved elements 9 (UCEs). We found that 194 (98.5%) of 197 non-exonic UCEs that are conserved across vertebrates align to the axolotl assembly. By comparison, 189 and 192 UCEs align to the Tibetan frog and Xenopus genomes, respectively, and 195 UCEs align to the coelacanth genome, indicating that the completeness of the axolotl genome assembly is comparable to or better than the two other amphibian genomes, which are smaller than 2.3 Gb 10 .
Salamanders serve as important tetrapod models for developmental, regeneration and evolutionary studies. An extensive molecular toolkit makes the Mexican axolotl (Ambystoma mexicanum) a key representative salamander for molecular investigations. Here we report the sequencing and assembly of the 32-gigabase-pair axolotl genome using an approach that combined long-read sequencing, optical mapping and development of a new genome assembler (MARVEL). We observed a size expansion of introns and intergenic regions, largely attributable to multiplication of long terminal repeat retroelements. We provide evidence that intron size in developmental genes is under constraint and that species-restricted genes may contribute to limb regeneration. The axolotl genome assembly does not contain the essential developmental gene Pax3. However, mutation of the axolotl Pax3 paralogue Pax7 resulted in an axolotl phenotype that was similar to those seen in Pax3 −/− and Pax7 −/− mutant mice. The axolotl genome provides a rich biological resource for developmental and evolutionary studies.
To further assess the completeness of the assembly, we generated a comprehensive gene catalogue by sequencing mRNA from 22 tissues (Supplementary Information section 3). Tissue-specific transcriptome assemblies and a composite assembly of all 1.5 billion transcript reads resulted in 180,649 transcript contigs (Supplementary Table 6) that contained 99% of the conserved core eukaryotic genes 11 and achieved the highest BUSCO score (http://busco.ezlab.org/) of an axolotl transcriptome reported to date (Supplementary Information section 3.4). More than 85% of the transcripts aligned to the genome along at least 95% of their length (Supplementary Information section 3.5), confirming the high completeness of the assembly. Furthermore, 71% of transcript contigs in which more than 95% of the sequence aligned with the genome were located on single scaffolds, demonstrating the high contiguity of the assembly. Using this comprehensive transcript set, we annotated a total of 23,251 protein-coding genes in the axolotl genome, a similar number to those found in other vertebrate genomes (Supplementary Information section 4.2).  Table 13) and included elements of more than 10 kb in length (Fig. 2c, Extended Data Fig. 1). Such long elements pose challenges for assembly, and indeed 97% of contigs ended in LTR elements. The number of substitutions to the repeat consensus, which is an estimate of the relative age of the LTR retroelement, indicates that the axolotl genome has undergone a long period of transposon activity followed by a recent and apparently continuing burst of expansion (Fig. 2d). This profile is consistent with previous small-scale characterizations of other salamander genomes 12 .

Expansion of long terminal repeat retroelements
The presence of many repeated elements contributes to a median intron size (22,759 bp) 13, 16 and 25 times that observed in human (1,750 bp), mouse (1,469 bp) and frog (906 bp), respectively (Fig. 3a, Supplementary Information section 4.3), a trend that was previously observed in five genes obtained from selective bacterial artificial chromosome sequencing of the axolotl genome 13 . Figure 3b shows a typical gene organization in axolotl compared to its human orthologue. Consistent with intron expansion, a distance comparison of pairs of highly conserved non-exonic elements shows that intergenic regions in the axolotl genome are 12 to 17 times larger than those in human, mouse and frog (Supplementary Information section 4.4).

HoxA cluster and intron size constraints
To examine gene cluster organization within this large genomic context, we focused on the HoxA locus, which has an important role in proximal-to-distal limb development and is reactivated during limb regeneration 14,15 . The entire HoxA locus is contained on a single contig (Fig. 3c), and the conserved neighbouring gene Evx1 is contained on the same 3.34-Mb scaffold. Compared to the orthologous human and frog clusters, the A. mexicanum HoxA cluster has a substantially increased repeat content and is 3.5 times larger, mostly owing to a 170-kb expansion between HoxA3 and HoxA4 (Fig. 3c). Notably, highly conserved non-exonic elements that putatively overlap cisregulatory elements are not interspersed in this 170-kb region, but remain in proximity to HoxA3 and HoxA4. The axolotl has a typical HoxA gene structure, with two coding exons separated by an intron. Notably, in contrast to the overall expansion of intron sizes, the intron sizes in the axolotl HoxA locus are very similar to those in other vertebrates, with the exception of AmHoxA3, which is also the longest of the HoxA genes in other tetrapods (Supplementary Table 17). Selected HoxC and HoxD genes examined in the red spotted newt exhibited similar properties 16 .

Article reSeArcH
On the basis of these observations, we examined the intron size distribution among a larger set of orthologous genes involved in developmental processes. While introns of non-developmental genes in axolotl show a median size expansion of 13-to 25-fold compared to human, mouse and frog, the expansion of introns of developmental genes is significantly lower (6-to 11-fold, P < 10 −11 ) (Fig. 3a, Supplementary Information section 4.3). In contrast to human, mouse and frog, introns of developmental genes in axolotl are shorter than introns of non-developmental genes. Furthermore, axolotl multi-exon genes that contain only short introns exhibit gene ontology enrichments related to developmental patterning that are not enriched in multi-exon genes with larger introns (Supplementary Table 16). These results suggest that intron size in developmental genes is under constraint in the axolotl, possibly because smaller gene sizes facilitate rapid transcription and thus upregulation of these genes in specific developmental contexts.

A reduced Pax-family complement
Next, we interrogated the genome for families of canonical developmental signalling molecules (Supplementary Information section 5). All three hedgehog paralogues as well as a full set of vertebrate Wnt genes were present (Extended Data Fig. 2a, b). However, we noted that certain members of the paired box family of transcription factors, which have diverse roles in tissue formation, were not found in the assembly. Consistent with the absence of Pax4 in amphibians and other vertebrate lineages 17 , the axolotl genome does not contain Pax4 but does contain Pax10. Notably, despite the presence of the Pax3 and Pax7 paralogues in all other known vertebrate lineages, we were able to identify Pax7 but not Pax3 in the axolotl genome assembly (Extended Data Fig. 2c). No Pax3 sequence was found in either the raw PacBio sequencing reads or the transcriptome. To confirm the loss of Pax3, we further examined the genomic region that would be syntenic for Pax3 for the presence of neighbouring genes and highly-conserved non-exonic elements (CNEs). The orthologues of genes surrounding mouse Pax3 (Sgpp2 and Epha4) were present in the A. mexicanum genome assembly; however, neither the Pax3 gene nor any of the Pax3associated CNEs were found (Fig. 3d). By contrast, several CNEs that overlap the Pax7 gene were identified in the assembly. Together, this evidence strongly suggests that Pax3 and several of its cis-regulatory elements are absent in the axolotl genome, probably owing to a deletion.

Axolotl Pax7 has similar functions to Pax3
To functionally assess the consequence of the absence of Pax3 in the axolotl, we used TALEN-and CRISPR-mediated gene editing 18 to mutate Pax7. In other vertebrates, Pax3 and Pax7 play key roles in muscle, neural tube and neural crest-derived tissue development 19 . Although these two genes share some common functions, deletion of Pax3 or Pax7 causes distinct phenotypes in mice [20][21][22] . We investigated whether frameshift deletions introduced into the AmPax7 gene would yield a comparable Pax7 phenotype, or whether AmPax7 may have taken on functions that are carried out by Pax3 in other vertebrates. Two different TALEN-mutant alleles (7-nt and 20-nt deletions) of AmPax7 were bred through two generations (Fig. 4a, Supplementary Information section 6). In the F2 generation, the developmental phenotype described below was observed in 83 out of 345 (24%) progeny from the Pax7 Δ20nt/+ intercrossing and 57 of 232 (24.6%) progeny from the Pax7 Δ7nt/+ intercrossing (Fig. 4b, Extended  Data Fig. 3). The phenotype was evident in homozygous mutants, as analysed by PCR and loss of protein (Supplementary Information sections 6.1, 6.3). This information, combined with the CRISPRmediated gene mutation results (Supplementary Information sections 6.2), shows that the homozygous Pax7 Δ20nt/Δ20nt and Pax7 Δ7nt/Δ7nt mutants represent recessive, complete or partial loss of Pax7 function.
The Pax7 Δ20nt/Δ20nt and strong F0 Pax7-CRISPR mutants exhibited a curved body, were unable to maintain an erect posture and exhibited a delay in growth compared to controls. Immunohistochemical analysis of trunk or tail cross-sections of early stage, 20-day-old Pax7 Δ20nt/Δ20nt or 17-day-old F0 Pax7-CRISPR axolotls showed normal muscle mass. However, at later stages, consistent with the mouse Pax7 deletion phenotype, tail and trunk muscles were greatly decreased (Fig. 4c, Extended Data Figs 4-6). Remarkably, the Pax7 mutant axolotls also completely lacked limb muscle (Figs 4d, Extended Data Fig. 7). In mice, Pax3, but not Pax7, is required for limb muscle formation [21][22][23] (Supplementary Table 18). These results demonstrate that AmPax7 has comparable functions to MmPax3 in the control of limb muscle genesis. In mice, Pax7 deletion affects craniofacial neural crest derivatives, including the facial bones 20 , whereas in zebrafish, pax7 mutants show loss of xanthophores and reduction of melanophores, but no loss of iridophores 24 . The AmPax7 mutants lacked a prefrontal bone, had a reduced number of melanophores and were deficient in xanthophores and iridophores except in the eyes (Fig. 4e-g, Extended Data Fig. 8). Pax3 deletion in mice is associated with neural tube closure defects 22,23 (Supplementary  Table 18). Similarly, Pax7 Δ20nt/Δ20nt -TALEN and Pax7-CRISPR axolotls displayed failed closure of the neural tube in the midbrain (Fig. 4h,  Extended Data Fig. 9). In summary, mutation of AmPax7 yields a combination of the Pax3-and Pax7-mutant phenotypes that are observed in other vertebrates ( Supplementary Information section 6). It will be interesting to understand how the regulation of Pax7 has changed in axolotl to enable the loss of Pax3, which is essential in other vertebrates.

Species-restricted genes in regeneration
Previous searches for mRNA and microRNA (miRNA) transcripts associated with limb regeneration relied on mapping to de novo transcriptome assemblies. We sought to re-examine these datasets using our newly acquired genomic data. Recent functional work has highlighted the role of diverged gene or protein function during regeneration [25][26][27] . Analysis of published tissue-enriched datasets 28 , combined with regeneration time courses 29,30 and our own transcriptional profiling of 22 tissues, identified five transcripts that are upregulated in the limb blastema (the mass of proliferating cells involved in regenerating the limb) with orthology limited to non-amniote vertebrates ( Supplementary Information section 7).
One of these protein sequences shows a weak similarity to tectorin, a basement membrane component normally found in the inner ear, consistent with studies that implicate extracellular matrix components with having an important role in limb regeneration 31,32 . Notably, another of these transcripts encodes a Ly6 family member in the urokinase type plasminogen activator surface receptor (uPAR) class. Previous studies had identified the salamander-specific Ly6 family member Prod1 as a key factor involved in salamander limb development and regeneration 25,33 . Our results suggest that Ly6 family members have a broader role in limb regeneration. Finally, we also investigated the role of non-coding RNAs by mapping a dataset of small RNA sequences expressed in the limb and limb blastema 34 to our genome assembly. This analysis classified 93 small RNAs as pre-miRNA sequences, of which 42 appear to be novel miRNAs ( Supplementary Information section 7.2). Taken together, these data point to a potential role in limb regeneration for several coding and non-coding sequences that have been lost or diverged rapidly 13  Article reSeArcH in amniotes. Future investigations of such sequences are likely to be a fruitful avenue for understanding the evolution of regeneration capabilities.

Discussion
We have generated a comprehensive whole-genome assembly for the salamander A. mexicanum, and analysis of this assembly has allowed us to draw conclusions about the structure of the expanded genome. Our data, together with data from plants and partial data from several other salamander species, show that LTR expansion is a major contributor to giant genome size across animals and plants 6,12,35 . Our assembly is sufficiently complete to reliably detect the absence of Pax3, which is present in fish and other amphibians. This analysis was confirmed using gene editing, which showed that AmPax7 has assumed functions that are carried out by Pax3 in other animals.
Functional analysis of axolotl development, physiology and regeneration is facilitated by the availability of tissue-and time-dependent gene expression profiles [28][29][30]36 . The axolotl genome provides a foundation for applying methods such as chromatin immunoprecipitation with sequencing (ChIP-seq) or assay for transposase-accessible chromatin using sequencing (ATAC-seq) to investigate the genomic basis of gene regulation during regeneration. Together with methods such as CRISPR-mediated gene editing, viral expression methods, transplantation and transgenesis, the axolotl is a powerful system for studying questions such as the evolutionary basis of its remarkable regeneration ability. Our approach of long-read sequencing, optical mapping and genome assembly using MARVEL also demonstrates that it is now feasible to assemble very large repeat-rich genomes.
Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper.

METhODS
No statistical methods were used to predetermine sample size. The experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment.
Axolotl genomic DNA was prepared from freshly isolated liver and spleen of an individual three year old adult d/d male using DNAzol followed by phenol/ chloroform extraction and ethanol precipitation.
A total of 50 size-selected SMRTbell libraries were prepared with a minimum fragment length cutoff between 10 kb and 20 kb. We sequenced medium and large insert libraries on the PacBio RSII instrument, making use of three different sequencing polymerases (P4, P5 and P6) and the corresponding sequencing chemistries (C2, C3 and C4). Movie times ranged from 180 min to 360 min with the majority of SMRT cells (1,414 of 1,933) at 240 min.
Sequences were assembled using the MARVEL assembler.
Optical mapping was performed using the Saphyr System (Bionano) based on NanoChannel array Technology. DNA was labelled with Nt.BspQI and Nb.BssSI enzymes in separate labelling reactions. Each enzyme reaction was run on the Saphyr System. 2.813 Tb of data were collected on three Saphyr Chips for Nt.BspQI and 2.0 Tb of data were collected on two Saphyr Chips for Nb.BssSI samples; single molecule N50 lengths were 240 kb and 184 kb, respectively. Each dataset was de novo assembled using Bionano Solve 2.1 software.
RNA was isolated from 22 tissue types using TRIzol or RNeasy reagents and sequenced using Illumina technology. The Trinity software package was used for transcriptome assembly. Code availability. The MARVEL assembler with documentation is available at https://github.com/schloi/MARVEL. Data availability. A browser of the axolotl genome is available at https://genome. axolotl-omics.org. The transcriptome assembly and the genome and transcriptome BLAST database can be accessed at https://www. axolotl-omics.org with no restrictions. The sequence data and both assemblies have been deposited in the NCBI BioProject database with accession numbers PRJNA378970 (genome data) and PRJNA378982 (transcriptome data). Both genome data and transcriptome data were deposited to the NCBI Nucleotide Database (nuccore) with accession numbers PGSH00000000 and GFZP00000000, respectively.

Data exclusions
Describe any data exclusions. We did not exclude data from the manuscript. In the final transcriptome assembly we did not include the oocyte dataset due to the presence of numerous short RNAs in that dataset that did not assemble with other contigs

Replication
Describe whether the experimental findings were reliably reproduced.
All attempts to replicate the mutant analysis were successful.

Randomization
Describe how samples/organisms/participants were allocated into experimental groups.
In mutant analysis, progeny were allocated as having a phenotype or no phenotype based on limb and body muscle mass. Then genotyping of both types of individuals was performed in the same experiment and the data examined for the genotype.

Blinding
Describe whether the investigators were blinded to group allocation during data collection and/or analysis.
The investigators were not blinded.
Note: all studies involving animals and/or human research participants must disclose whether blinding and randomization were used. For manuscripts utilizing custom algorithms or software that are central to the paper but not yet described in the published literature, software must be made available to editors and reviewers upon request. We strongly encourage code deposition in a community repository (e.g. GitHub). Nature Methods guidance for providing algorithms and software for publication provides further information on this topic.

Materials and reagents
Policy information about availability of materials

Materials availability
Indicate whether there are restrictions on availability of unique materials or if these materials are only available for distribution by a for-profit company.
All unique materials are openly available.

Antibodies
Describe the antibodies used and how they were validated for use in the system under study (i.e. assay and species).  Up to juvenile stages, it is not possible to determine the sex of the axolotls by morphology and no cytochemical or molecular assay is available. The phenotype is very likely independent of the gender.