Chromosome-scale genome assembly reveals how repeat elements shape non-coding RNA landscapes active during newt limb regeneration

Abstract

Newts have large genomes harboring many repeat elements. How these elements shape the genome and relate to newts' unique regeneration ability remains unknown. We present here the chromosome-scale assembly of the 20.3 Gb genome of the Iberian ribbed newt, Pleurodeles waltl, with a hitherto unprecedented contiguity and completeness among giant genomes. Utilizing this assembly, we demonstrate conserved synteny as well as genetic rearrangements, such as in the major histocompatibility complex locus. We provide evidence suggesting that intronic repeat elements drive newt-specific circular RNA (circRNA) biogenesis and show their regeneration-specific expression. We also present a comprehensive in-depth annotation and chromosomal mapping of microRNAs, highlighting genomic expansion profiles as well as a distinct regulatory pattern in the regenerating limb. These data reveal links between repeat elements, non-coding RNAs, and adult regeneration and provide key resources for addressing developmental, regenerative, and evolutionary principles.

Keywords: PacBio; amphibian; blastema; circular RNA; giant genome; microRNA; non-coding RNA; regeneration; salamander; transposable element.

Graphical abstract

Figure 1

Sequencing and chromosome-level assembly of the P. waltl genome (A) Schematic representation of the sequencing and assembly strategy. (B) P. waltl genome assembly features (top) and BUSCO assessment (bottom). (C) Contig N(X) plot showing what percentage of each assembled genome (X) is contained within pieces at least N(X) bp in size. Shown are contig statistics from P. waltl (this study and Elewa et al.¹¹), Protopterus annectens, and Ambystoma mexicanum. (D) Hi-C interaction heatmap of contact data for scaffolded genome. Individual scaffolds are delineated. Denser areas of red signal off diagonal represent interactions between the arms of the same chromosome. (E) Gene completeness based on BUSCO single-copy Vertebrata orthologs (n = 3,354). Scores are based on annotations for P. waltl, P. annectens, and A. mexicanum.

Figure 2

Comparative synteny of the P. waltl genome (A) Oxford plot depicting the locations of chordate linkage groups (CLGs). A total of 14,062 best-reciprocal-hit orthologs were identified between the P. waltl chromosomes and 17 ancestral CLGs. The dense rectangular blocks of dots represent units of deeply conserved synteny. Fisher’s exact test (FET) was used to calculate the significance of interactions between each scaffold and linkage group. Dots are colored by CLG, with solid dots indicating the FET p ≤ 0.05 and translucent dots depicting the FET p > 0.05. (B) Ribbon plot showing syntenic blocks of the newt, axolotl and gar genomes. Best-reciprocal-hit orthologs are connected by ribbons between the P. waltl, the A. mexicanum, and the L. oculatus genomes and colored based on their identification as proteins from the 17 CLGs as in (A). (C and D) Regions of the newt, axolotl, and gar genomes identified pairwise as containing co-localized blocks of genes with the CLGs between axolotl and newt (C) and gar and newt (D). The size of each circle corresponds to the number of genes in each identified syntenic block. Circles are colored by CLG. Individual numbers are also available in Table S2. Regions of the P. waltl genome identified as containing syntenic blocks of genes from the chordate linkage groups (CLG), related to Figure 2, Table S3. Regions of the A. mexicanum genome identified as containing syntenic blocks of genes from the chordate linkage groups (CLG), related to Figure 2, Table S4. Regions of the L. oculatus genome identified as containing syntenic blocks of genes from the chordate linkage groups (CLG), related to Figure 2, Table S5. Regions of the P. waltl genome identified as containing syntenic blocks of genes from the A. mexicanum genome, related to Figure 2, Table S6. Regions of the P. waltl genome identified as containing syntenic blocks of genes from the L. oculatus genome, related to Figure 2. (E) Microsynteny analysis of the Fgf5 locus between gar, lungfish, frog, newt, and axolotl. The newt locus has experienced various disruptions, including a 58 Mb gap between Prdm8 and Cfap299 as well as an inversion of order between Paqr3 and Prdm8. Arrows represent the relative position and direction of genes on chromosomes, retrogenes are marked by white rectangles, and filled rectangles imply a different chromosomal position.

Figure 3

Genome expansion, repeat element composition, and regeneration-associated expression of hATs and domesticated hATs in P. waltl (A) Bar plots showing the contribution of the main repeat element classes to each of the indicated genomes. Note the abundance of DNA repeat elements in P. waltl. (B) Heatmap representing percentage contribution of repeat element families to the genome of each indicated species. Families are grouped by repeat element type. (C) Relative contributions of hAT families to the total repeat element repertoire of the P. waltl genome, expressed as a percentage of all repeat element based on sequence length. (D) Expression of indicated hAT elements in PacBio Iso-seq-derived transcriptomes from P. waltl limb blastema, brain, and spleen. (E) Differential expression of hATs (left) as well as domesticated hATs, Zbed1, Zbed4, Zbed6, and Zbed8 (right), during limb regeneration in P. waltl, based on normalized and centered RNA-seq counts for the indicated conditions. Only genes whose differential expression is significantly altered (between 0 and 3 dpa or 0 and 7 dpa) are depicted. Color key represents Z-score values.

Figure 4

Characterization of the circular RNA transcriptome in P. waltl (A) Venn diagram of intersection between different circRNA prediction programs. A total of 9,799 circRNA candidates were detected by all three programs. (B) Longer introns flank predicted circRNAs compared to introns flanking randomly generated backspliced junctions. The mean length of up- and downstream flanking introns for predictions was 49,805 nt and for random permutations was 29,619 nt. U Int, upstream flanking introns; D Int, downstream flanking introns; ∗∗∗∗p < 0.0001 (Welch’s t test). (C) Bar plot of “hotspot” genes, which host 15 or more circRNA isoforms. (D) Schematic of circRNA formation promoted by inverted-repeat elements in flanking introns (top). Table of repeat elements that frequently overlap with down- and upstream introns flanking P. waltl circRNAs (bottom). (E) Principal-component analysis plot based on circRNA expression across tissue groups. Tissue groups include adult tissue (brain, eyes, heart, liver, and lung), adult forelimb (n = 2), adult limb stump (D0 limb stump tissue), developing larvae (larvae at limb bud stage, late embryo stages 22 + 25), and adult regenerating limb (n = 2 of 3 dpa and n = 3 of 7 dpa). dpa, days post-amputation. (F) Volcano plot of differentially expressed circRNAs in adult regenerating limb (3 and 7 dpa) compared to developing larvae (limb bud stage and s22 + s25). dpa, days post-amputation.

Figure 5

Annotation and genomic organization of miRNAs in P. waltl (A) Pipeline for annotation of miRNAs. (B) Heatmap of normalized mature miRNA read counts across tissues for all miRDeep2 candidates. Forty-three ESCC-miRNAs are highlighted. (C) Pie chart of the genomic localization for all 855 miRNAs and the 43 ESCC-miRNAs. (D) Bar plot representing the percentage of miRNAs that are embedded in or located within 100 bp of a repeat element. (E) Distribution of all the P. waltl miRNAs at the chromosome level. All 855 miRNA genomic loci are indicated on the left and the 43 ESCC-miRNAs are indicated on the right. ESCC, embryonic stem cell-specific cell-cycle regulating.

Figure 6

Conserved miRNA paralogs and multiple alignment analyses of P. waltl ESCC-miRNAs (A) Heatmap of conserved miRNA paralogs across species of various genome sizes. Color scale corresponds to the number of paralogs (log2) in each organism, and arrows indicate miR-427/430/302 and P. waltl. (B) Unrooted phylogenetic tree of the ESCC-miRNA precursors. Each spoke represents one ESCC-miRNA precursor, where greater distance between spokes indicates greater sequence diversity. (C) MAFFT alignment of all 43 ESCC-miRNAs containing the conserved AAGUGC seed site on either the 5′ or the 3′ arm. Darker blue indicates higher conservation of bases. MAFFT, multiple alignment using fast Fourier transform.