Recherche

... Article https://doi.org/10.1038/s41467-023-43012-9 Conserved chromatin and repetitive patterns reveal slow genome evolution in frogs Received: 28 October 2021 Accepted: 27 October 2023 1234567890():,; 1234567890():,; Check for updates Jessen V. Bredeson 1,2,19, Austin B. Mudd1,19, Soa Medina-Ruiz1,19, Therese Mitros1, Owen Kabnick Smith 3, Kelly E. Miller1, Jessica B. Lyons 1, Sanjit S. Batra4, Joseph Park1, Kodiak C. Berkoff 1, Christopher Plott 5, Jane Grimwood 5, Jeremy Schmutz 5, Guadalupe Aguirre-Figueroa3, Mustafa K. Khokha 6, Maura Lane6, Isabelle Philipp1, Mara Laslo 7, James Hanken 7, Gwenneg Kerdivel8, Nicolas Buisine8, Laurent M. Sachs 8, Daniel R. Buchholz9, Taejoon Kwon 10,11, Heidi Smith-Parker 12, Marcos Gridi-Papp 13, Michael J. Ryan12, Robert D. Denton14, John H. Malone 14, John B. Wallingford15, Aaron F. Straight 3, Rebecca Heald 1, Dirk Hockemeyer 1,16,17, Richard M. Harland 1 & Daniel S. Rokhsar 1,2,16,17,18 Frogs are an ecologically diverse and phylogenetically ancient group of anuran amphibians that include important vertebrate cell and developmental model systems, notably the genus Xenopus. Here we report a high-quality reference genome sequence for the western clawed frog, Xenopus tropicalis, along with draft chromosome-scale sequences of three distantly related emerging model frog species, Eleutherodactylus coqui, Engystomops pustulosus, and Hymenochirus boettgeri. Frog chromosomes have remained remarkably stable since the Mesozoic Era, with limited Robertsonian (i.e., arm-preserving) translocations and end-to-end fusions found among the smaller chromosomes. Conservation of synteny includes conservation of centromere locations, marked by centromeric tandem repeats associated with Cenp-a binding surrounded by pericentromeric LINE/L1 elements. This work explores the structure of chromosomes across frogs, using a dense meiotic linkage map for X. tropicalis and chromatin conformation capture (Hi-C) data for all species. Abundant satellite repeats occupy the unusually long (~20 megabase) terminal regions of each chromosome that coincide with high rates of recombination. Both embryonic and differentiated cells show reproducible associations of centromeric chromatin and of telomeres, reecting a Rabl-like conguration. Our comparative analyses reveal 13 conserved ancestral anuran chromosomes from which contemporary frog genomes were constructed. Amphibians are widely used models in developmental and cell biology15, and their importance extends to the elds of infectious disease, ecology, pharmacology, environmental health, and biological diversity610. While the principal model systems belong to the genus A full list of afliations appears at the end of the paper. Nature Communications | (2024)15:579 Xenopus (notably the diploid western clawed frog X. tropicalis and the paleo-allotetraploid African clawed frog X. laevis), other amphibian models have increasingly been introduced due to their diverse developmental, cell biological, physiological, and behavioral adaptations1121. e-mail: dsrokhsar@gmail.com 1 Article While genome evolution has been extensively studied in mammals22 and birds23,24, the relative lack of phylogenetically diverse chromosome-scale frog genomes has limited the study of genome evolution in anuran amphibians. Here, we report a high-quality assembly for X. tropicalis and three new chromosome-scale genome assemblies for the Puerto Rican coqu (Eleutherodactylus coqui), a direct-developing frog without a tadpole stage16,19, the tngara frog (Engystomops pustulosus), which is a model for vocalization and mate choice15,18,20, and the Zaire dwarf clawed frog (Hymenochirus boettgeri), which has an unusually small embryo, is a model for regulation of cell and body sizes, and a source of potent host-defense peptides with therapeutic potential13,17,21. Genome assemblies are essential resources for further work to exploit the experimental possibilities of these diverse animals. The new high-quality X. tropicalis genome upgrades previous draft assemblies25,26 and our new genomes complement draft chromosome-scale sequences for the African clawed frog27 (Xenopus laevis), the African bullfrog28 (Pyxicephalus adspersus), the Leishan moustache toad29 (Leptobrachium leishanense), the Ailao moustache toad30 (Leptobrachium [Vibrissaphora] ailaonicum), and Asiatic toad31 (Bufo gargarizans), as well as scaffold- and contig-scale assemblies for other species32. The rapidly increasing number of chromosome-scale genome assemblies makes anurans ripe for comparative genomic and evolutionary analysis. Chromosome number variation among frogs is limited3335. Based on cytological36,37 and sequence comparisons25,27,33,38,39 most frogs have n ~1012 pairs of chromosomes. A recent meiotic map of the yellowbellied toad Bombina variegata showed that its twelve chromosomes are simply related to the ten chromosomes of X. tropicalis40. The stability of the frog karyotype contrasts with the more dramatic variation seen across mammals22,37,41,42, which as a group is considerably younger than frogs. The constancy of the frog karyotype parallels the static karyotypes of birds23,43, although birds typically have nearly three times more chromosomes than frogs, including numerous microchromosomes (among frogs, only the basal Ascaphus44 has microchromosomes). Despite the stable frog chromosome number, however, fusions, ssions, and other interchromosomal rearrangements do occur, and we can use comparisons among chromosomescale genome sequences to (1) infer the ancestral chromosomal elements, (2) determine the rearrangements that have occurred during frog phylogeny, and (3) characterize the patterns of chromosomal change among frogs. These ndings of conserved synteny among frogs are consistent with prior demonstrations of conservation between Xenopus tropicalis with other tetrapods, including human and chicken25,45. Since frog karyotypes are so highly conserved, X. tropicalis can be used as a model for studying chromosome structure40, chromatin interaction, and recombination for the entire clade. Features that can be illuminated at the sequence level include the structure and organization of centromeres and the nature of the unusually long subtelomeres relative to mammals (frog subtelomeres are ~20 megabases, compared with the mammalian subtelomeres that are typically shorter than a megabase). The extended subtelomeres of frogs form interacting chromatin structures in interphase nuclei that reect threedimensional intra-chromosome and inter-chromosome subtelomeric contacts, which are consistent with a Rabl-like conguration. As in other animals, subtelomeres of frogs have an elevated GC content and recombination rate. Here we show that the unusually high enrichment of recombination in the subtelomeres likely reects similar structural and functional properties in other vertebrates, though the quality of the assembly reveals that the length of subtelomeres, expansion of microsatellite repeat sequences by unequal crossing over, and high recombination rates are considerably greater in frogs than in mammals. A strong correlation between recombination rate and microsatellite sequences suggests that unequal crossing over during meiotic recombination is implicated in the expansion of satellites in the Nature Communications | (2024)15:579 https://doi.org/10.1038/s41467-023-43012-9 subtelomeres. We use Cenp-a binding at satellites to conrm centromere identity and extend the predictive power of the repeat structures to centromeres of other frogs. We address the unusually high recombination rate in subtelomeric regions, correlating with the landscape of base composition and transposons. Over the 200 million years (My) of evolution that we address here, centromeres have generally been stable, but the few karyotypic changes reveal the predominant Robertsonian translocations at centromeric regions; we also document the slow degeneration that occurs to inactivated centromeres and fused telomeres, changes that are obscured in animals with rapidly evolving karyotypes. Results and discussion High-quality chromosome-scale genome assembly for X. tropicalis To study the structure and organization of Xenopus tropicalis chromosomes and facilitate comparisons with other frog genomes, we assembled a high-quality chromosomal reference genome sequence (Supplementary Data 1, Supplementary Fig. 1, and Supplementary Notes 1 and 2) by integrating data from multiple sequencing technologies, including Single-Molecule Real-Time long reads (SMRT sequencing; Pacic Biosciences), linked-read sets (10x Genomics), short-read shotgun sequencing, in vivo chromatin conformation capture, and meiotic mapping, combined with previously generated dideoxy shotgun sequence. New sequences were generated from 17thgeneration individuals from the same inbred Nigerian line that was used in the original Sanger shotgun sequencing45. The new reference assembly, version 10 (v10), spans 1448.4 megabases (Mb) and is substantially more complete than the previous (v9) sequence25, assigning 219.2 Mb more sequence to chromosomes (Supplementary Table 1). The v10 assembly is also far more contiguous, with half of the sequence contained in 32 contigs longer than 14.6 Mb (in comparison, this N50-length was. 71.0 kilobases [kb] in v9). The assembly captures 99.6% of known coding sequences (Supplementary Table 2 and Supplementary Note 2). We found that the fragmented quality of earlier assemblies was due, in part, to the fact that 68.3 Mb (4.71%) of the genome was not sampled by the 8 redundant Sanger dideoxy whole-genome shotgun dataset45 (Supplementary Fig. 2ac and Supplementary Note 2). These missing sequences are apparently due to non-uniformities in shotgun cloning and/or sequencing (Supplementary Fig. 2df). Previously absent sequences are distributed across 140.5k blocks of mean size 485.7 basepairs (bp) (longest 50.0 kb) on the new reference assembly, are enriched for sequences with high GC content (Supplementary Fig. 2g), and capture an additional 6774 protein-coding exons from among 4718 CDS sequences (Supplementary Fig. 2d, e). The enhanced contiguity of v10 is accounted for by the relatively uniform coverage of PacBio long-read sequences along the genome, as expected from other studies4649. Most remaining gaps are in highly repetitive and satellite-rich centromeres and subtelomeric regions (see below) (Supplementary Fig. 2a). Additional chromosome-scale frog genomes To assess the evolution of chromosome structure across a diverse set of frogs, we generated chromosome-scale genome assemblies for three new emerging model species, including the Zaire dwarf clawed frog Hymenochirus boettgeri (a member of the family Pipidae along with Xenopus spp.), and two neobratrachians: the Puerto Rican coqu Eleutherodactylus coqui (family Eleutherodactylidae) and the tngara frog Engystomops pustulosus (family Leptodactylidae). These chromosome-scale draft genomes were primarily assembled from short-read datasets and chromatin conformation capture (Hi-C) data (Supplementary Data 1, Supplementary Table 3, and Supplementary Note 3). To further expand the scope of our comparisons, we also updated the assemblies of two recently published frog genomes: the 2 Article African bullfrog Pyxicephalus adspersus28, from the neobatrachian family Pyxicephalidae, and the Ailao moustache toad Leptobrachium (Vibrissaphora) ailaonicum29, from the family Megophryidae (Supplementary Fig. 3 and Supplementary Note 3). These species span the pipanuran clade, which comprises all extant frogs except for a small number of phylogenetically basal taxa, such as Bombina40 and Ascaphus50. The chromosome numbers of the new assemblies agree with previously described karyotypes for E. coqui51 (2n = 26) and E. pustulosus52 (2n = 22). The literature for H. boettgeri, however, is more equivocal, with reports53,54 of 2n = 2024. The n = 9 chromosomes of our H. boettgeri assembly are consistent with our chromosome spreads (Supplementary Fig. 3a). The karyotype variability in the published literature and discrepancy with the karyotypes of our H. boettgeri samples may be the result of cryptic sub-populations within this species or segregating chromosome polymorphisms. Protein-coding gene set for X. tropicalis The improved X. tropicalis genome encodes an estimated 25,016 protein-coding genes (Supplementary Table 4), which we predicted by taking advantage of 8580 full-length-insert X. tropicalis cDNAs from the Mammalian Gene Collection55 (MGC), 1.27 million Sangersequenced expressed sequence tags45 (ESTs), and 334.5 gigabases (Gb) of RNA-seq data from an aggregate of 16 conditions and tissues56,57 (Supplementary Data 1 and Supplementary Note 2). The predicted gene set is a notable improvement on previous annotations, both in completeness and in full-length gene-level accuracy, due in part to the more complete and contiguous assembly (Supplementary Fig. 1, Supplementary Table 2, and Supplementary Note 2). In particular, singlemolecule long reads lled gaps in the previous X. tropicalis genome assemblies that likely arose from cloning biases in the Sanger sequencing process, encompassing exons embedded in highly repetitive sequences (Supplementary Fig. 2). A measure of this completeness and the utility of the X. tropicalis genome is provided by comparing its gene set with those of vertebrate model systems with reference-quality genomes, including chicken58, zebrash59, mouse60, and human61,62 (Supplementary Fig. 4ac). Notably, despite the closer phylogenetic relationship between birds and mammals, X. tropicalis shares more orthologous gene families (and mutual best hits) with human than does chicken, possibly because of the loss of genomic segments in the bird lineage23,63 and/or residual incompleteness of the chicken reference sequence, due to the absence of several microchromosomes58. For example, of 13,008 vertebrate gene families with representation from at least four of the vertebrate reference species, only 341 are missing from X. tropicalis versus 1110 from chicken (Supplementary Fig. 4a). The current X. tropicalis genome assembly also resolves gene order and completeness of gene structures in the long subtelomeres that were missed in previous assemblies due to their highly repetitive nature (Supplementary Fig. 2). Protein-coding gene sets for additional frogs We annotated the new genomes of E. coqui, E. pustulosus, H. boettgeri, and P. adspersus using transcriptome data from these species (Supplementary Data 1) and peptide homology with X. tropicalis (Supplementary Tables 5 and 6). To include mustache toad in our cross-frog comparisons, we adopted the published annotation from ref. 29 (Supplementary Note 3). We found 14,412 orthologous groups across the ve genera with OrthoVenn264, including genes found in at least four of the ve frog genera represented (Supplementary Fig. 4d). As expected, due to its reference-quality genome and well-studied transcriptome, only 72 of these clusters were not represented in X. tropicalis (and only 42 clusters from gene families present in six or more members among a larger set of seven frog species, see Supplementary Fig. 4e); the additional frog genomes each had between 575 and 712 of these genes missing (or mis-clustered), suggesting better than 95% Nature Communications | (2024)15:579 https://doi.org/10.1038/s41467-023-43012-9 completeness in the other species. For analyses of synteny, we further restricted our attention to 7292 one-to-one gene orthologs that were present on chromosomes (as opposed to unlinked scaffolds) in the core genomes X. tropicalis, H. boettgeri, E. coqui, E. pustulosus, and P. adspersus. The total branch length in the pipanuran tree shown in Fig. 1 (including both X. laevis subgenomes) is 2.58 substitutions per fourfold synonymous site. Repetitive landscape Centromeric and telomeric tandem repeats play a critical role in the stability of chromosome structure65. Nonetheless, other kinds of repeats also play a role in the preservation of these important chromosome landmarks66, 67. The new X. tropicalis v10 assembly captures sequences from centromeres and distal subtelomeres that were fragmented in the previous assemblies25,45. The percentage of the genome covered by transposable elements is slightly higher than previously reported45 (36.82% vs. 34%) (Supplementary Table 7). Insertional bias in the pericentromeric regions is observed for specic families of long interspersed elements (LINEs), including the relatively young Chicken Repeat 1 (CR1)68 (3.14% of the genome) and the ancient L1 (1.06%) (Fig. 2 and Supplementary Fig. 5). The X. tropicalis v10 assembly captures signicantly more tandem repeats in the distal subtelomeric portions of the genome relative to earlier assemblies. An exhaustive search for tandem repeats using Tandem Repeats Finder69 determined that 10.67% of the chromosomes are covered by tandem arrays consisting of 5 or more monomeric units greater than 10 bp. Many tandem repeat footprints lie in the gaps of previous assemblies25,45 (Supplementary Fig. 2). Our new hybrid genome assembly closed many gaps containing centromeric and subtelomeric tandem repeats, and captured numerous subtelomeric genes (Supplementary Fig. 2). The overall repeat landscape derived from the X. tropicalis assembly is mirrored in the other frog assemblies, with similar centromeric repeats, and lengthy subtelomeres, as discussed below. Genetic variation The inbred X. tropicalis reference genotype was nominally derived from 17 generations of brother-sister mating, starting with two Nigerian founders. In the absence of selection, this process should lead to an increasingly homozygous genome due to increasing identity by descent of the two reference haplotypes, with residual heterozygosity conned to short blocks totaling a fraction ~1.17 (0.809)t of the genetic map70, or 3.2% after t = 17 generations of full-sib mating. In contrast, we observe that 11.7% of the genome (125.12 cM out of a total of 1070.16 cM) exhibits residual heterozygosity (Supplementary Fig. 6). While this excess could be explained by balancing selection due to recessive lethals, a more mundane possibility is that some non-fullsib mating occurred during the inbreeding process. Errors early in the inbreeding process would be consistent with the unexpectedly high heterozygosity (~44%) observed in two 13th-generation members of the lineage (Supplementary Fig. 6), which far exceeds the 7.4% theoretical expectation from repeated full-sib mating. The approximately fourfold further reduction from these individuals to our 17thgeneration reference, however, is consistent with theoretical expectations in the absence of selection. Residual blocks of heterozygosity after inbreeding reect distinct founder haplotypes. Within these blocks, we observe 3.0 singlenucleotide variants per kilobase, which serves as an estimate of the heterozygosity of the wild Nigerian population. To begin to develop a catalog of segregating variation in X. tropicalis, we also shotgunsequenced pools of frogs from the Nigerian and Ivory Coast B populations, which are the two main sources of experimental animals. These two populations have been previously analyzed using SSLP markers71. From our light pool shotgun analysis, we identied a total of 6,546,379 SNPs, including 2,482,703 variants in the Nigerian pool and 4,661,928 3 Article https://doi.org/10.1038/s41467-023-43012-9 2n = 26 2n = 26 2n = 22 2n = 26 2n = 26 2n = 26 2n = 26 2n = 26 2n = 18 2n = 20 2n = 20 2n = 20 2n = 18 2n = 18 Mya 2n = 18 Fig. 1 | Phylogenetic tree and gene ortholog alignment. The phylogenetic tree of the seven analyzed species, calculated from fourfold degenerate sites and divergence time condence intervals, drawn with FigTree (commit 901211e, https:// github.com/rambaut/gtree): Xenopus tropicalis, X. laevis, and Hymenochirus boettgeri (Pipoidea: Pipidae); Leptobrachium (Vibrissaphora) ailaonicum (Pelobatoidea: Megaphrynidae); Engystomops pustulosus (Neobatrachia [Hyloidea]: Leptodactylidae), Eleutherodactylus coqui (Neobatrachia [Hyloidea]: Euleutherodactylidae); and Pyxicephalus adspersus (Neobatrachia [Ranoidea]: Pyxicephalidae). The ancestral karyotype is labeled at each node on the tree. Black circles with white text refer to chromosome changes summarized in Table 1. The alignment plot was generated with JCVI using the 7292 described chromosome one- to-one gene orthologs from OrthoVenn2, followed by manual ltering of single stray orthologs. The Hi-C-derived centromere position is represented with a black circle on each chromosome. Ancestral chromosomes (A to M) are labeled at the top of the alignment based on the corresponding region in P. adspersus. The alignments for each ancestral chromosome are colored uniquely, with those upstream and downstream of the X. tropicalis centromeric satellite repeat colored in dark and light shades of the ancestral chromosome color. Chromosomes labeled with asterisks are shown reverse complemented relative to their orientations in the genome assembly. Mya millions of years ago, n the haploid chromosome number. Source data are provided as a Source Data le. in the Ivory Coast B pool, with 598,252 shared by both pools, suggesting differentiation between populations (Supplementary Fig. 6 and Supplementary Note 2). At least some of these pipanuran elements have a deeper ancestry within amphibians. For example, the chromosomes of the discoglossid frog Bombina variegata (n = 12), an outgroup to the pipanurans, show considerable conservation of synteny with X. tropicalis based on linkage mapping40. Compared with the pipanuran ancestral elements described here, the nine B. variegata chromosomes 2, 3, 4, 5, 6, 8, 9, 10, and 12 correspond to nine pipanuran elements A, B, C, F, G, H, I, E, and J, respectively, extending these syntenic elements to the last common ancestor of Bombina+pipanurans (which does not have a common name). The remaining three B. variegata chromosomes 1, 7, and 11 are combinations of the remaining four pipanuran elements D, K, L, and M. Similarly, the genome of the axolotl, Ambystoma mexicanum, a member of the order Caudata (salamanders and newts) and ~292 million years divergent from pipanurans74, also conserves multiple syntenic units with pipanurans (Supplementary Fig. 7i). For example, axolotl chromosomes 4, 6, 7, and 14 are in near 1:1 correspondence with pipanuran elements F, A, B, and K, respectively, although small pieces of F and A can be found on axolotl 10, and parts of B can be found on axolotl 9 and 13. Other axolotl chromosomes are fusions of parts of two or more pipanuran elements. For example, axolotl chromosome 5 is a fusion of a portion of J with most of G; the remainder of G is fused with a portion of L on the q arm of axolotl chromosome 2. Further Conserved synteny and ancestral chromosomes Comparison of the chromosomal positions of orthologs across seven frog genomes reveals extensive conservation of synteny and collinearity (Fig. 1 and Supplementary Fig. 7ag). We identied 13 conserved pipanuran syntenic units that we denote A through M (Methods and Supplementary Note 4). Each unit likely represents an ancestral pipanuran chromosome, an observation consistent with the 2n = 26 ancestral karyotype inferred from cytogenetic comparisons across frogs36,72. Over 95% (6952 of 7292) of chromosomal one-to-one gene orthologs are maintained in the same unit across the ve frog species, attesting to the stability of these chromosomal elements (Fig. 1). The conservation of gene content per element is comparable to the 95% ortholog maintenance in the Muller elements in Drosophila spp73. Despite an over twofold difference in total genome size across the sampled genomes, each ancestral pipanuran element accounts for a nearly constant proportion of the total genome size, gene count, and repeat count in each species, implying uniform expansions and contractions during the history of the clade (Supplementary Fig. 7h). Nature Communications | (2024)15:579 4 https://doi.org/10.1038/s41467-023-43012-9 Repeat enrichment 150 75 5.4 2.7 Rec. Rate kb/Mb Article Chr1 Chr2 Chr3 Chr4 Chr5 Chr6 Chr7 Chr8 Chr9 Chr10 Fig. 2 | Density of pericentromeric and subtelomeric repeats in Xenopus tropicalis. Pericentromeric (red) and subtelomeric (purple) regions were used to obtain enriched repeats, excluding chromosomes with short p-arms (chromosomes 3, 8, and 10). Pericentromeric repeats (yellow) correspond to selected subsets of non-LTR retrotransposons (CR1, L1, and Penelope), LTR retrotransposons (Ty3), and DNA transposons (PiggyBac and Harbinger). Subtelomere- enriched repeats (blue) correspond mainly to satellite repeats and LTR retrotransposons (Ty3, Ngaro). Densities of each repeat type plotted as kb/Mb. Chromosomes are centered by the position of centromeric tandem repeats (black dots). Rates of recombination (Rec. rate) in cM/Mb are shown as solid black lines. Tick marks indicate 10 Mb blocks (Supplementary Fig. 5). kb kilobases, Mb megabases, cM centiMorgans. Source data are provided as a Source Data le. comparisons are needed to determine which of these rearrangements occurred on the axolotl vs. the stem pipanuran lineage. Genomes from the superfamilies Leiopelmatoidea and Alytoidea, which diverged prior to the radiation of pipanurans, will also be informative. Chromosomal conserved synteny across pipanuran frogs is comparable to that observed in birds, which have evolved by limited intra-chromosomal rearrangement from an n = 40 ancestor43, mostly involving fusion of microchromosomes, as we nd here for pipanurans (see below). The relative stasis of frog and bird chromosomes is in contrast to the variable karyotypes of mammals, which was rst noted by Bush et al.37 and is now extensively documented at the level of chromosomal painting22 and genome sequence42. The reasons for these different modes of evolution remain unclear but are likely related to the difculty in xing partial-arm chromosomal rearrangements in large historically panmictic populations due to reduced fertility in translocation heterozygotes, as rst noted by Wright75. Partial-arm rearrangements, as observed in mammals, can become xed in populations that are dynamically subdivided by local extinction and colonization, which allows the reduced fertility of translocation heterozygotes to be overcome by genetic drift76. Robertsonian or centric translocations involving breaks and joins near centromeres account for several of the rare rearrangements (Figs. 1 and 3b). For example, element G clearly experienced centric ssion in the E. coqui lineage. Conversely, I and M underwent centric fusion in the E. pustulosus lineage. E. coqui has experienced the most intense rearrangement, including Robertsonian ssions of A and G, a Robertsonian fusion of I/K, and a signicant series of Robertsonian rearrangements involving B, E, F, and H that resulted in Bprox/H, Bdist/ Fdist, and E/Fprox (Table 1 and Supplementary Table 8). (Mechanistically, these ssions and fusions likely occur by translocations; see ref. 77 for a discussion.) Elements I and H form the two arms of a submetacentric chromosome in pipids (Fig. 3a), and therefore the pipid ancestor, but are found as either independent acrocentric chromosomes (e.g., in P. adspersus and L. ailaonicum) or as arms of Table 1 | Organization and conservation of the 13 ancestral chromosomes of pipanuran genomes Phylogenetic position Structural event (1) Stem pipid lineage J + K JK D. + E. D.E Chromosome evolution Block rearrangements of the 13 ancestral elements dominate the evolutionary dynamics of pipanuran karyotypes (Table 1 and Fig. 1). While element C has remained intact as a single chromosome across the group (except for internal inversions), all of the other elements have experienced translocations during pipanuran evolution. During these translocations, the elements have remained intact except for the breakage of elements A and M by reciprocal partial-arm exchange observed in P. adspersus chromosomes 3 and 6. To trace the evolutionary history of centromeres shown in Fig. 1, we inferred their positions using Hi-C contact map patterns, as in X. tropicalis (where centromeres were also conrmed by analysis of Cenp-a binding as described below). In general, the pericentromeres of other pipanurans were characterized by the same repetitive element families found in Xenopus, further corroborating their identication. Overall, we found broad pericentromeric conservation among the species analyzed (Figs. 1 and 3a). Nature Communications | (2024)15:579 I + H I H (Rob. fusion) (2) P. adspersus lineage after divergence from R. temporaria (3) E. pustulosus lineage after divergence from E. coqui (4) E. coqui lineage after divergence from E. pustulosus A + M A1.m1 + m2.A2 M + I M.I (Rob) K + D K.D (Possible end-end) G1 G2 G1 + G2 (Rob. ssion) A1 A2 A1 + A2 (Rob. ssion) I + K I K (Rob. fusion + inversion) E + F1F2 + B1B2 + H EF1 + F2B2 + B1H (5) H. boettgeri lineage after divergence from Xenopus M + JK MJK (6) X. laevis progenitor lineage after divergence from X. tropicalis L + M LM Rob Robertsonian. Middle-dots (i.e., ) represent centromeres. Periods (i.e., .) represent translocation breakpoints. 5 Article https://doi.org/10.1038/s41467-023-43012-9 a Hbo 7 0 Xtr 7 0 10 Xla 7L 0 10 Epu 1 0 Xtr 1 0 Eco 7 0 20 40 Lai 10 0 20 40 Xtr 4 0 20 40 60 20 30 20 80 40 30 100 50 40 50 120 60 140 70 60 80 70 160 180 200 90 100 110 80 90 100 220 120 110 240 130 120 b 20 40 60 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 60 80 100 120 Eco 3 80 100 120 0 140 20 160 180 200 40 60 80 100 120 140 160 40 60 80 100 120 140 160 c Pad 9 10 60 60 20 40 Lai 8 0 80 100 120 30 20 40 50 60 0 70 80 20 90 100 Pad 8 0 110 120 20 130 40 140 150 60 Fig. 3 | Subtelomeric repeats highlight regions of chromosome fusion. Examples of (a) conserved structure and pericentromere maintenance of H. boettgeri (Hbo), X. tropicalis (Xtr), and X. laevis (Xla) chromosomes; b a Robertsonian translocation in the lineage leading to E. coqui (Eco), shown compared with E. pustulosus (Epu) and X. tropicalis; and c an end-to-end fusion that occurred in the lineage giving rise to X. tropicalis and subsequent pericentromere loss, shown compared with L. ailaonicum (Lai) and P. adspersus (Pad). The analyzed species were visualized with a custom script, alignment_plots.py (v1.0, https://github.com/ abmudd/Assembly). For each plot, the Hi-C inference-based centromeric regions are depicted with black stars, the X. tropicalis centromeric satellite repeat from tandem repeat analysis with a red star (on X. tropicalis chromosomes 7 and 1 (a, b), the stars overlap), the density of L1 repeats per chromosome with gold densities, and the runs of collinearity containing at least one kilobase of aligned sequence between the species with connecting black lines. kb kilobases, Mb megabases. Source data are provided as a Source Data le. (sub)metacentrics formed by centric fusion with other elements (Supplementary Table 8). We also observed end-to-end fusions78 of (sub)metacentric chromosomes, for example, the joining of D with K in E. pustulosus, and with element E in the common ancestor of pipids (Hymenochirus and Xenopus) (Figs. 1 and 3c). Since bicentric chromosomes are not stably propagated through mitosis, one of the two ancestral centromeres brought together by end-to-end fusion must be lost or inactivated, as shown in Fig. 3c for the ancient DE fusion in pipids. We note that the D centromere persists in both end-to-end fusions involving D, suggesting that centromeres derived from different ancestral elements may be differentially susceptible to silencing, although with only two examples this could have happened by chance. Using the pericentromeric and subtelomeric repeats landscape as a proxy, we found several examples of end-to-end chromosome fusions in which residual subtelomeric signals are preserved near the presumptive junctions (Fig. 3 and Supplementary Fig. 8). These include the end-to-end fusion of X. tropicalis-like chromosomes 9 and 10 (elements L and M) to produce the X. laevis chromosome 9_10 progenitor that is found in both the L and S subgenomes of this allotetraploid27. These X. laevis chromosomes display evidence of decaying subtelomeric signatures in the region surrounding the ancestral LM fusion (Fig. 1 and Supplementary Fig. 8a, b). Similarly, enrichment of subtelomerically-associated repeats is observed in H. boettgeri chromosome 8_10 (Supplementary Fig. 8ce) near the junction between the portions of the chromosome with M and J/K ancestry (the J/K fusion occurred near the base of pipids). In both cases, the centromere from element M (i.e., the centromere in X. tropicalis chromosome 9) is maintained after fusion. The inversion of the p-arm from chromosome 8S also has evidence of decaying sequence but the median is less than the median Jukes-Cantor (JC) distance at the chromosome 9_10 fusion, suggesting that the fusion preceded the inversion. Rate of karyotype change Nature Communications | (2024)15:579 The long-range and, in most cases, chromosome-scale collinearity (Supplementary Fig. 7 and Supplementary Table 9) among the frog species we examined, despite a combined branch length of 1.05 billion years (Supplementary Tables 10 and 11), parallels the conserved synteny observed in birds79 and reptiles80, but differs from the substantial chromosome variation found in mammals22,41. Maintenance of collinear blocks may reect an intrinsically slow rate of rearrangement in frogs, perhaps a consequence of large regions devoid of recombination, or selection favoring retention of specic gene order and chromosome structure related to chromosomal functions. We inferred 8 fusions, 2 ssions, one pairwise, and one four-way reciprocal fusion; counting the last as a composite of three pairwise rearrangements yields a total of 17 translocations (excluding smaller intra-chromosome rearrangements) corresponding to an average rate of one karyotype change every 62 million years (Fig. 1 and Table 1). This rate is similar to the rate of one chromosome number change every 70 to 90 million years as previously proposed for frogs and some mammals33,37 but still slower than karyotype change rates for most mammals81 and many reptiles82. Of course, our rate calculation is based on only seven species, and the rate may vary depending on the species analyzed. Some frog taxa, such as Eleutherodactylus spp. (2n = 1632) and Pristimantis spp51. (2n = 2238), have experienced higher rates of karyotype change. On the other hand, other lineages, such as those leading to Leptobrachium ailaonicum, L. leishanense14, and Rana temporaria83, have had no detectable inter-chromosome exchange over the past 205 million years (Fig. 1). Nonetheless, this analysis of chromosome variation across the frog lineage is consistent with an overall slow rate of karyotype evolution84. Considering rearrangement rate variation across taxa, we can ask whether any of the individual branches show an unusually high or low number of translocations relative to the overall pipanuran rate. The absolute karyotype stasis of L. ailaonicum over ~200 My is only 6 Article marginally slower than the pipanuran average (two-sided test, P = 0.04 under a simple Poisson model of 1 change every 62 My, before familywise correction for testing of multiple lineages). Conversely, the E. coqui lineage has experienced six translocations during a time interval in which only one rearrangement would be expected. This is a signicant enrichment relative to the Poisson model (P = 1 103) and is the only branch on which the constant rate hypothesis is rejected. Notably, Euleutherodactylus is the most karyotypically variable frog genus, suggesting possible ongoing karyotypic instability84,85. Regarding chromosome stability, our collection only includes one example in which a chromosome arm is disrupted by translocation; all other changes are either Robertsonian (involving breaks near a centromere) or end-to-end (near a telomere). This observation allows us to reject (P < 4 104) a simple random break model, under which we would expect ~12.3 chromosome arms to be broken across our phylogeny (Supplementary Note 4). This suggests that centromeric and telomeric regions are more prone to breakage, and/or breaks within chromosome arms are selected against. The latter model is consistent with a reduced probability of xation of reciprocal (partial-arm) translocations due to selection against reduced fertility in heterozygotes75, which can be overcome by genetic drift under some conditions76. Centromeres, satellites, and pericentromeric repeats The stasis of Xenopus chromosomes relative to other frogs (see above) allows us to examine the repetitive landscape of chromosomes that are not frequently rearranged by translocation and may be approaching a structural equilibrium. Vertebrate centromeres are typically characterized by tandem families of centromeric satellites (e.g., the alpha satellites of humans) that bind to the centromeric histone H3 protein, Cenp-a, a centromere-specic variant of histone H365,86. Cenp-a binding satellites have been described in X. laevis87, and here we nd distantly related X. tropicalis satellite sequences that also co-precipitate with Cenp-a. Thus, chromatin immunoprecipitation and sequencing (ChIP-seq) shows that Cenp-a binding coincides with the predictions of centromere positions derived from chromatin conformation analysis and repetitive content (Supplementary Figs. 5ac and 9ac and Supplementary Tables 12 and 13). Importantly, this concordance supports the prediction of centromere position for other species that we infer below. The Cenp-a-bound sequences are arrays of 205-bp monomers that share a mean sequence identity greater than 95% at the nucleotide level, with a specic segment of the repeating unit showing the greatest variability (Supplementary Fig. 9d, e). The X. tropicalis centromere sequence is different from centromeric-associated repeats found in X. laevis87,88, suggesting the sequences evolve rapidly after speciation but are maintained across chromosomes within the species. All pericentromeric regions of (sub)metacentric X. tropicalis chromosomes are enriched in retrotransposable repetitive elements (15 Mb regions shown in Fig. 2). In other vertebrate species and Drosophila, retrotransposable elements from the pericentromeric regions are involved in the recruitment of constitutive heterochromatin components89,90. Among the pericentromerically-enriched repeats we identied specic families belonging to LTR retrotransposons (Ty3), non-LTR retrotransposons (CR1, Penelope, and L1), and DNA transposable elements (PIF-Harbinger and piggyBac families) (Fig. 2 and Supplementary Fig. 5). CR1 (CR1-2_XT) is the most prevalent and among the youngest of all pericentromeric retrotransposons (mean JukesCantor (JC) distance to consensus of 0.05). In contrast, L1 and Penelope types have a mean JC greater than 0.4 (Supplementary Fig. 5). The age of the repeats, indirectly measured by the JC distance, suggests that pericentromeric retrotransposons have experienced different bursts of activity and tendency to insert near the centromere. Expression of active retrotransposons and random insertion can compromise chromosome stability, and because silencing of these is crucial, genomes develop mechanisms to rapidly silence them. Such Nature Communications | (2024)15:579 https://doi.org/10.1038/s41467-023-43012-9 insertions may be positively selected, and therefore amplied, to establish pericentromeric heterochromatin, but may be counterselected when they insert in gene-rich chromosome arms. Recombination and extended subtelomeres With chromosome sequences in hand, we studied the distribution of recombination along X. tropicalis chromosomes using a previously generated Nigerian-Ivory Coast F2 cross25 (Supplementary Note 5 and Supplementary Data 2). Half of the observed recombination is concentrated in only 160 Mb (11.0% of the genome) and 90% of the observed recombination occurs in 540 Mb (37.3%). In contrast, the extended central regions of each chromosome are cold, with recombination rates below 0.5 cM/Mb and that are often indistinguishable from zero in our data (Supplementary Fig. 10a, b and Supplementary Table 14). Strikingly, we nd that (sex-averaged) recombination is concentrated within just 30 Mb of the ends of each chromosome and occurs only rarely elsewhere (Supplementary Fig. 10a). The regions of the subtelomeres experiencing high recombination are nearly sixfold longer than in non-amphibian genomes91,92. The rates of recombination in Xenopus subtelomeres were not previously determined, since the repeat-rich subtelomeres were absent from earlier assemblies, and markers present in those regions showed insufcient linkage to be incorporated into linkage maps25. Elevated rates of recombination near telomeres and long central regions of low recombination have been observed in the macrochromosomes of diverse tetrapods, including birds92,93, snakes94, and mammals9597. This pattern appears to be independent of the involvement of the chromatin modier PRDM9 in dening recombination hotspots98 since dogs lack PRDM9 but show the same pattern, with elevated recombination in promoter regions and around CpG islands96. Conversely, snakes possess the prdm9 gene but also show hotspots of recombination concentrated in promoters and functional regions94. Since amphibians lack the prdm9 gene99, we further analyzed the genomic features that colocalized in subtelomeric regions prone to recombination. To assess sequence features associated with enriched recombination, we focused on the extended subtelomeres, dened as the terminal 30 Mb of all (sub)metacentric chromosomes and the terminal 30 Mb excluding the 15 Mb surrounding the pericentromeric regions of acrocentric chromosomes (3, 8, and 10) (Fig. 2). The median recombination rate in the extended subtelomeres (1.72 cM/Mb) is over tenfold higher than the median rate observed in the rest of the chromosome arms (0.14 cM/Mb) (two-sample KolmogorovSmirnov test, two-sided, Hochberg-corrected P = 5.2 10321) (Supplementary Fig. 10c and Supplementary Note 5). The recombination rate in the 5-Mb region surrounding the centromeric tandem repeats is even lower (0.01 cM/Mb). Since constitutive heterochromatin in pericentromeric regions is known to repress recombination, this observation is expected (reviewed in refs. 100,101). However, the centromeres of acrocentric chromosomes lie within 30 Mb of telomeres and preclude the presence of extended subtelomere-associated repeats (Fig. 2 and Supplementary Fig. 11). We examined the relationship between rates of recombination against repetitive elements and sequence motifs associated with recombination hotspots in other vertebrate species (Supplementary Fig. 12a and Supplementary Table 14). Similar to chicken and zebra nch, recombination is the highest in subtelomeres and positively correlates with GC content92,93,102, which is consistent with GC-biased gene conversion83,103,104 in recombinogenic regions (median GC = 42.5% in the 74 Mb in which half of the recombination occurs) vs. the nonrecombinogenic centers of chromosomes (median 38.8%). As in zebra nch (Supplementary Fig. 13), recombination in X. tropicalis is strongly correlated with satellite repeats (Pearsons correlation, r = 0.68, R2 = 0.457). The high density of satellite repeats (Supplementary Table 15) in highly recombinogenic subtelomeric regions suggests that 7 Article https://doi.org/10.1038/s41467-023-43012-9 unequal crossing over during meiotic recombination mediates tandem repeat expansions105,106. Notably, in the extended subtelomeric regions tandem repeats are enriched in specic tetrameric sequences (TGGG, AGGG, and ACAG) compared to non-tandem repeats (Supplementary Fig. 12b). In contrast, centromeric tandem repeats are completely devoid of these short sequences. Some of the tandem arrays enriched in the terminal 30 Mb of all chromosomes derive from portions of transposable elements, such as SINE/tRNA-V, LINE/CR1, DNA/Kolobok-2 (Supplementary Fig. 11 and Supplementary Table 16). For example, the minisatellite expansion that arose from the family of SINE/tRNA-V present in the pipid lineage107 amplied a 52-bp portion of the 3UTR-tail from the SINE/ tRNA-V element in Xenopus tropicalis and other frog species (Supplementary Table 17). Although intact SINE/tRNA-V elements are distributed throughout the genome, the minisatellite fragment is only expanded in subtelomeric SINE/tRNA-Vs, suggesting that recombination in subtelomeres has driven minisatellite expansion (Supplementary Figs. 11 and 14). Interestingly, although the satellite expansions are similar in X. laevis and X. tropicalis, they differ in other frogs, suggesting that different satellite expansions can occur repeatedly during the maintenance of the long subtelomeric regions (see below). We hypothesize that the high rate of recombination in the extended subtelomeres of frog chromosomes drives tandem repeat expansion through illegitimate homologous recombination and, in the process, increases GC content (Supplementary Fig. 14d, e). Unfortunately, it is difcult to resolve cause and effect with observational data, To further rene our understanding of chromosome structure in X. tropicalis, we studied chromatin conformation capture (Hi-C) data from nucleated blood cells. These experiments link short reads representing sequences in close three-dimensional proximity108. Figure 4 shows mapped Hi-C read pairs for chromosomes 1 and 2, with different minimum mapping quality thresholds above and below the diagonal (Supplementary Fig. 1e and Supplementary Note 5). We consistently observe a wing of intra-chromosome contacts transverse to the main diagonal, which (1) intersects the main diagonal near the cytogenetically dened Cenp-a-binding centromere, and (2) indicates contacts between p and q-arms (Supplementary Figs. 1e and 15). These observations imply that interphase chromosomes are folded at their centromeres, with contacts between distal arms. We also observe enriched inter-chromosome contacts among centromeres and among chromosome arms along a centromere-to-telomere axis, suggesting that chromosomes are organized in a polarized arrangement in the nucleus (Supplementary Figs. 9a and 15 and Supplementary Table 18). Notably, the correlation between centromere position and the observed intra-chromosome folding and inter-chromosome contacts at centromeres allows us to use Hi-C analysis and principal Fig. 4 | Organization of X. tropicalis chromosomes into Rabl-like conguration and distinct nuclear territories. a Hi-C contact matrices for chromosomes 1 and 2 (lower-left and upper-right gold boxes, respectively) showing features of the threedimensional chromatin architecture within X. tropicalis blood cell nuclei. Blue pixels represent chromatin contacts between XY pairs of 500 kb genomic loci, with intensity proportional to contact frequency. Hi-C read pairs are mapped stringently (MQ 30) above the diagonal and permissively (MQ 0) below the diagonal. The characteristic A/B-compartment (checkerboard) and Rabl-like (angel wing) interarm contact patterns within each chromosome are evident. Above the diagonal, an increased frequency of interchromosomal chromatin contacts is observed between pericentromeres (connected by dotted lines) and between chromosome arms (Supplementary Tables 18, 19, and 21), suggesting a centromere-clustered organization of chromosomes in a Rabl-like conguration. Below the diagonal, high-intensity pixels near the ends of chromosomes not present above the diagonal suggest a telomere-proximal spatial bias in the distributions of similar genomic repeats. See Supplementary Fig. 1e for a plot showing all chromosomes. b Chromosome territories within the nucleus. Yellow, white, and blue colors indicate the normalized relative enrichment, parity, and depletion of chromatin contacts between non-homologous chromosomes (Supplementary Tables 21 and 22). For example, chromosome 1 exhibits higher relative contact frequencies with all chromosomes except chromosomes 7, 9, and 10, which are generally depleted of contacts except among themselves (MQ 30; 2 (81, n = 24,987,749) = 3,049,787; Hochberg-corrected P < 4.46 10308; Relative range: 0.827741.16834). Note, due to the inbred nature of the Nigerian strain, contacts could not be partitioned by haplotype, and so the results reported here represent chromosomal averages. c Schematic representation of chromosome territories from (b). The size of each chromosome number is approximately proportional to the number of enriched interactions. Darker and lighter colors indicate chromosomes nearer and more distant to the reader, respectively. Mb megabases, MQ mapping quality. Source data are provided as a Source Data le. Nature Communications | (2024)15:579 and we cannot rule out the alternative hypothesis that meiotic recombination is promoted by preferential DNA breakage at short sequence motifs (Supplementary Fig. 12b), which is then repaired by homologous recombination. Chromatin conformation correlates with cytogenetic features 8 Article component analysis (PCA) of intra- and inter-chromosome contacts109 to infer the likely centromeric positions based purely on Hi-C data in frogs whose cytogenetics are less well-studied (see below). Taken together, these intra- and inter-chromosome contacts in Xenopus blood cells are consistent with a Rabl-like (Type-I110) chromosome conguration111, 112. Such associations among centromeres and among telomeres, rst observed in salamander embryos111, have been observed in other animals110,113117, fungi110,118,119, and plants109,110,120122. Outside of mammals, Rabl-like contacts have been observed in a wide diversity of taxa. Hoencamp et al.110. surveyed 24 plant and animal species using Hi-C and observed Rabl-like patterns in 14 (58.3%) of them. Out of seven vertebrates sampled, however, only Xenopus laevis broblasts showed a Rabl-like pattern. We note that Hi-C patterns can depend on cell type, cell cycle stage, and developmental time; and while Rabl-like Hi-C patterns are often absent from tissue samples used in mammalian genome sequencing projects, they have been observed in studies of mouse and human cell lines (Supplementary Note 5). In X. tropicalis, this conguration is understood to be a relict structure from the previous mitosis123,124 in which the chromosomes have become elongated and telomeres clustered on the inner nuclear periphery. Dernburg and colleagues125 reasoned that the Rabl conguration observed in Drosophila embryonic nuclei126,127 is a result of anaphase chromosome movement and, due to their rapidly dividing nature, such chromosomes are unable to relax into a diffused chromatin state. Consistent with this, we nd that Rabl-like chromosomal interarm contacts in early frog development (NF stages 823) appear more tightly constrained (mean SEM: sum of squared distances [SSD] 1.384 0.066, centromere-to-telomere-polar interarm contact enrichment [CTP] 2.492 0.179) in these rapidly dividing cells. Notably, more specialized (liver and brain) X. tropicalis adult tissues, except for blood cell nuclei (SSD 1.465, CTP 1.813), show less chromosomal interarm constraint (mean SEM: SSD 5.233 1.258, CTP 1.362 0.153) (Supplementary Fig. 16, Supplementary Table 19, and Supplementary Note 5). Although it is possible that some amount of Hi-C signal may be due to residual incompleteness in the assembly and concomitant mismapping of reads to repeat sequences, these observations are robust to quality ltering, even when using singlecopy sequences. Furthermore, such contacts are similarly weak in sperm cells16 (SSD 6.285, CTP 1.056), a control that argues strongly against sequence mismapping artifacts (Supplementary Note 5). As noted above, the presence and strength of Rabl-like congurations vary depending on the tissue, cell type, and developmental time. Such variability highlights the need to sample a broader diversity of tissues and time points to characterize completely the Rabl-like chromosome structures in X. tropicalis. Chromatin compartments Chromatin contacts in human108,128,129, mouse129, chicken130 and other phylogenetically diverse species131133 often show a characteristic checkerboard pattern that is superimposed on the predominant neardiagonal signal. This pattern implies an alternating A/B-compartment structure with enriched intra-compartment contacts within chromosomes (Fig. 5a), which has been linked with G-banding in humans134. X. tropicalis also exhibits an A/B-compartment pattern, which emerges as alternating gene-rich (A) and gene-poor (B) regions (median 19.99 genes/Mb and 9.99 genes/Mb, respectively) (Fig. 5b). Despite their twofold difference in gene content, A and B-compartment lengths are comparable, with approximately exponential distributions (Supplementary Fig. 17). The arithmetic mean sizes are A = 1.32 Mb, B = 1.48 Mb; the corresponding geometric means (i.e., the exponential of the arithmetic mean of logarithms of lengths) are somewhat shorter (A = 0.807 Mb, B = 0.946 Mb). A/B compartments are also differentiated by repetitive content129, with A-compartment domains showing slight enrichment (1.211.44-fold) in DNA transposons of the DNA/Kolobok-T2, DNA/hAT-Charlie, and Mariner-Tc1 families. Nature Communications | (2024)15:579 https://doi.org/10.1038/s41467-023-43012-9 B-compartment domains had signicantly higher enrichment for DNA transposons (DNA/hAT-Ac, Mar-Tigger) and retrotransposons (Ty3/ metaviridae and CR1), among other repeats (1.122.11-fold) (Fig. 5c, Supplementary Table 20). The association between repeats overrepresented in A and B compartments is also captured in one of the principal components obtained from the repeat densities of all chromosomes (Supplementary Note 5); we detect a modest negative correlation (Pearsons r = 0.44) between A/B compartments and the third principal component obtained from the repeat density matrix (Supplementary Fig. 5b). The association between chromatin condensation and repeat type could be due to a preference for certain transposable elements to insert in specic chromatin contexts, or chromatin condensation to be controlled, in part, by transposable element content, or a combination of these factors. However, we were unable to nd any correlation of A/B compartments with the G-banding of condensed chromosomes in X. tropicalis135,136. Higher-order chromatin interactions Chromatin conformation contacts also provide clues to the organization of chromosomes within the nucleus. We observe non-random (2 (81, n = 24,987,749) = 3,049,787; Hochberg-corrected P < 4.46 10308) associations between chromosomes in blood cell nuclei (Fig. 4b and Supplementary Tables 21 and 22): (a) chromosome 1 is enriched for contacts with chromosomes 28 (mean 1.05 enrichment), and depleted of contacts with 9 and 10 (mean 0.89); (b) among themselves, chromosomes 28 show differential contact enrichment or depletion; and (c) chromosomes 9 and 10 are enriched (1.17) for contacts with one another, but are depleted of contacts with all other chromosomes. These observations suggest the presence of distinct chromosome territories111,137139, where chromosomes 28 are localized more proximal toand arrayed aroundchromosome 1, with chromosomes 9 and 10 relatively sequestered from chromosome 1 (Fig. 4c). The contact enrichment between chromosomes 9 and 10 is particularly notable because these short chromosomes (91.2 and 52.4 Mb, respectively) have become fused in the X. laevis lineage140, which might have been enabled by their persistent nuclear proximity141143. Between chromosomes, p-p and q-q arm interactions exhibit a small but signicant enrichment (1.059 enrichment; 2 (1, n = 24,786,496) = 17,037; Hochberg-corrected P < 4.46 10308) over p-q arm contacts. This is a general feature of (sub)metacentric chromosomes observed in other frog genomes (Supplementary Table 21), except E. coqui (0.928 enrichment; 2 (1, n = 6,850,547) = 3,914; Hochberg-corrected P < 4.46 10308), the chromosomes of which appear predominantly acrocentric or telocentric. Finally, the p-arms of chromosomes 3, 4, 8, and 9 are enriched for contacts with both p and q-arms of chromosome 10, with the acrocentric chromosomes 3 and 8 showing the strongest relative enrichment and a slight preference between p-arms. The q-arms of chromosomes 3 and 8, however, exhibit a slight enrichment for contacts with the larger (sub)metacentric chromosomes 1, 2, 4, and 5. Taken together, these observations suggest possible colocalization of the p and q-arms of chromosomes 3 and 8 in X. tropicalis blood cell nuclei. Future impacts Anuran amphibians play a central role in biology, not simply as a globally distributed animal group, but also as key subjects for research in areas that range from ecology and evolution to cell and developmental biology. The genomic resources generated here will thus provide important tools for further studies. Given the crucial role of X. tropicalis for genomic analysis of development and regeneration144,145, the improvements to our understanding of its genome reported here will provide a more nely-grained view of biomedically important genetic and epigenetic mechanisms. This new genome is also important from the standpoint of evolutionary genomics, as comparisons between the genomes of X. tropicalis and X. laevis shed light on the 9 Article c 210 Compartment Genes / Mb 200 190 180 170 160 100 25 A 120 d 110 B 1 2 3 4 5 6 7 8 9 10 100 90 70 60 400 200 0.10 0.05 0.00 0.05 0.10 DNA/hAT-Ac DNA DNA/TcMAr-Tigger LTR/Ty3-Metaviridae LINE1/CR1 Ty3-Metaviridae DNA_transposon hAT DIRS 21 10 0 210 200 190 180 170 160 150 140 130 120 110 100 80 90 70 60 50 40 30 20 0 10 0.15 e Eigenvector (HiC) 0.15 1 Eigenvector (PC3: Repeat density) 210 200 190 180 170 160 150 140 130 120 110 80 90 100 70 60 50 40 30 20 0 10 10 Transposable_Element 20 LTR_Retrotransposon 30 Mariner/Tc1 DNA/Kolobok-T2 40 DNA/TcMar-Tc1 0 50 0 Repeat class 600 Repeats / Mb Position (Mb) 0 130 80 Eigenvector 50 50 140 b A B 75 0 150 Chromosome 100 DNA/hAT-Charlie a https://doi.org/10.1038/s41467-023-43012-9 0 Position (Mb) 20 40 60 80 100 120 140 160 180 200 Position (Mb) Fig. 5 | A/B-compartment structure and gene/repeat densities. a Correlation matrix of intra-chromosomal Hi-C contact densities between all pairs of nonoverlapping 250 kb loci on chromosome 1. Yellow and blue pixels indicate correlation and anti-correlation, respectively, and reveal which genomic loci occupy the same or different chromatin compartment. Black pixels indicate weak/no correlation. b The rst principal component (PC) vector revealing the compartment structure along chromosome 1, obtained by singular value decomposition of the correlation matrix in panel a. Yellow (positive) and blue (negative) loadings indicate regions of chromosome 1 partitioned into A and B compartments, respectively. c Gene density (genes per megabase) distributions in A (yellow) vs. B (blue) compartments genome-wide and per chromosome. Sample sizes and signicance statistics provided in Supplementary Table 20. d Repeat classes signicantly enriched by density (repeats per megabase) in A (yellow) vs. B (blue) compartments. Sample sizes and signicance statistics provided in Supplementary Table 20. Each boxplot summarizes the combined (A + B) density distribution (Y-axis) per class (X axis); lower and upper bounds of each box (black) delimit the rst and third quartiles, respectively, and whiskers extend to 1.5 times the interquartile range, while the median per class is represented as a lled white circle. e The PC3 loadings (purple line) from the repeat density matrix inversely correlate with alternating A/B-compartment loadings (green) for chromosome 1. See Supplementary Fig. 5b for all chromosomes. Purple rectangles plotted on the X axis denote subtelomeric regions, the red rectangle spans the pericentromere, and the black point marks the median centromere-associated tandem repeat position. Mb megabases. Source data are provided as a Source Data le. consequences of genome duplication145. The new genome described here for H. boettgeri, another pipid frog, is also signicant in this regard, as it enables an interesting comparison of Xenopus genomes to that of a closely related outgroup. Moreover, the genomes of E. coqui and E. pustulosus provide a foundation for future studies of the evolution of ontogenies and their underlying developmental mechanisms, as E. coqui is a direct-developing frog with no tadpole stage16 and E. pustulosus, a foam-nesting frog, is a model for studying mating calls and female mate choice18. In addition to their interesting life histories, both frogs display distinct patterns of gastrulation146,147. Finally, recent work has demonstrated the efcacy of genetic or genomic analysis for understanding the impact of chytrid fungus on various amphibian species148. A deeper and broader understanding of amphibian genomes will be useful in the context of the global decline of amphibian populations149,150. Note added in proof: The recent nding of tetraploid dwarf clawed frogs from the Congo suggests that the diploid Hymenochirus we studied may distinct from H. boettgeri151. Xenopus tropicalis genomic DNA extraction and sequencing Methods This study complies with the ethical standards set forth by the Institutional Animal Care and Use Committee (IACUC) protocols at the University of California Berkeley, Yale University, University of Cincinnati, and the University of the Pacic. The IACUC and associated facilities are subject to review and oversight by NIHs Ofce of Lab Animal Welfare. Nature Communications | (2024)15:579 High molecular weight DNA was extracted from the blood of an F17 Xenopus tropicalis Nigerian strain female25. Paired-end (PE) Illumina whole-genome shotgun (WGS) libraries were constructed by the QB3 Functional Genomics Laboratory (FGL) using a KAPA HyperPrep Kit and sequenced on an Illumina HiSeq 2500 as 2 250 bp reads by the Vincent J. Coates Genomics Sequencing Lab (VCGSL) at the University of California, Berkeley (UCB). Single-Molecule Real-Time (SMRT) continuous long-read (CLR) sequencing was performed at the HudsonAlpha Institute for Biotechnology (HAIB) on Pacic Biosciences (PacBio) RSII machines with P6-C4 chemistry (Supplementary Data 1). Chromium Genome linked-read (10x Genomics) sequencing was carried out by HAIB on an Illumina HiSeq X Ten. Hi-C libraries were constructed by Dovetail Genomics LLC. See Supplementary Note 1 for more detailed extraction and sequencing methods. Xenopus tropicalis genome assembly and annotation Chromium linked-read (10x Genomics) data were assembled with Supernova152 (v1.1.5). This assembly was used to seed the assembly of PacBio CLR data using DBG2OLC153 (commit 1f7e752). An independent PacBio-only assembly was constructed with Canu154 (v1.6-132-gf9284f8). These two assemblies were combined, or metassembled, using MUMmer155 (v3.23) and quickmerge156 (commit e4ea490) (Supplementary Fig. 1a). Residual haplotypic redundancy was identied and removed (Supplementary Fig. 1b). The non-redundant metassembly 10 Article was scaffolded with Sanger paired-ends and BAC-ends45 using SSPACE157 (v3.0) and Hi-C using 3D-DNA117,158,159 (commit 2796c3b), then manually curated in Juicebox160,161 (v1.9.0). The assembly was polished with Arrow162 (smrtlink v6.0.0.47841), Pilon163 (v1.23), and then FreeBayes164 (v1.1.0-54-g49413aa) with ILEC (map4cns commit dd89f52, https:// bitbucket.org/rokhsar-lab/map4cns). The genome was annotated with the DOE-Joint Genome Institute (JGI) Integrated Gene Call (IGC) pipeline165 (v5.0) using transcript assemblies (TAs) generated with Trinity166,167 (v2.5.1) from multiple developmental stages and tissues (Supplementary Data 1). RepeatModeler168 (v1.0.11) was run on all frog species. The frog and ancestral repeat libraries from RepBase169 (v23.12) were combined with the repeat consensuses identied by RepeatModeler. The merged repeat library was used to annotate repeats of all frogs with RepeatMasker170 (v4.0.7). See Supplementary Note 2 for more detailed assembly and annotation methods. Hymenochirus boettgeri metaphase chromosome spread H. boettgeri were obtained from Albany Aquarium (Albany, CA). Stage 26 tadpoles (n = 10) were incubated at room temperature in 0.01% colchicine and 1 MMR for 46 h. After removing the yolky ventral portion of the tadpoles, the remaining dorsal portions were pooled together in deionized water and allowed to stand for 20 min. The dorsal portions were transferred to 0.2 mL of 60% acetic acid in deionized water and allowed to stand for 5 min. The tissue was then pipetted onto a positively charged microscope slide, and excess acetic acid was blotted away. To atten the tissue and promote chromosome spreading, the slide was covered with a coverslip, and a lead brick was placed on top of it for 5 min. The slide and coverslip were then placed on dry ice for 5 min. The coverslip was removed from the frozen slide, and the slide was stained with 0.1 mg/mL Hoechst Stain solution for 5 min. A fresh coverslip was then mounted on the slide using VectaShield, and the edges were sealed with nail polish. Chromosomes in metaphase spreads (Supplementary Fig. 3a) were imaged on an Olympus BX51 Fluorescence Microscope run with Metamorph (v7.0) software using a 60 oil objective. Chromosome number was counted in 75 separate metaphase spreads. Genome and transcriptome sequencing of ve pipanurans Illumina PE 10x Genomics Chromium linked-read whole-genome libraries for E. pustulosus (from liver), E. coqui (from blood), and H. boettgeri (from liver) were sequenced on an HiSeq X at HAIB. PacBio SMRT Sequel I CLR data were generated at UC Davis DNA Technologies and Expression Analysis Core for each of E. pustulosus and H. boettgeri from liver samples. In addition, two Illumina TruSeq PE libraries (from kidney) and two Nextera mate-pair libraries (from liver) for E. coqui were prepared. Hi-C libraries were prepared for H. boettgeri, E. pustulosus, and E. coqui using the DovetailTM Hi-C Kit for Illumina (Beta v0.3 Short manual) following the Animal Tissue Samples protocol, then sequenced on a HiSeq 4000 at the VCGSL or a NextSeq at Dovetail Genomics. Illumina TruSeq Stranded mRNA Library Prep Kit (cat# RS-122-2101 and RS-122-2102) libraries were prepared from E. pustulosus stages 45 and 56 whole tadpoles (gut excluded) and various adult tissues dissected from frogs maintained at the University of the Pacic. Brain (n = 3), dorsal skin (n = 2), eggs (n = 2), eye (n = 2), heart (n = 2), intestine (n = 2), larynx (n = 3), liver (n = 2), lung (n = 2), and ventral skin (n = 2) samples were washed twice with PBS, homogenized in TRIzol Reagent, and centrifuged, followed by ash freezing of the supernatant. RNA was isolated following the TRIzol Reagent User Guide (Pub. No. MAN0001271 Rev. A.0) protocol. In addition, H. boettgeri eggs were homogenized in TRIzol Reagent and processed according to the manufacturers instructions. RNA was then isolated using the QIAGEN RNeasy Mini Kit (cat# 74104). An Illumina mRNA library was prepared using the Takara PrepX RNA-Seq for Illumina Library Kit (cat# 640097) by the QB3 FGL at UCB. All libraries were sequenced at the VCGSL on an Nature Communications | (2024)15:579 https://doi.org/10.1038/s41467-023-43012-9 HiSeq 4000 as 2 151 bp reads. See Supplementary Note 3 for additional details about DNA/RNA extractions and library preparations, and Supplementary Data 1 for a complete list of DNA/RNA sequencing data generated for E. coqui, E. pustulosus, and H. boettgeri. Assembly and annotation of ve pipanuran genomes E. pustulosus and H. boettgeri contigs were assembled with Supernova152 (v2.0.1). E. coqui contigs were assembled with Meraculous171,172 (v2.2.4) and residual haplotypic redundancy was removed using a custom script (align_pipeline.sh v1.0, https://github. com/abmudd/Assembly) before scaffolding with SSPACE157 (v3.0). E. pustulosus and H. boettgeri contigs were ordered and oriented using MUMmer155 (v3.23) alignments to PBEC-polished (map4cns commit dd89f52, https://bitbucket.org/rokhsar-lab/map4cns) DBG2OLC153 (commit 1f7e752) hybrid contigs (Supplementary Note 3). All three assemblies were scaffolded further with linked reads and Scaff10X (v2.1, https://sourceforge.net/projects/phusion2/les/scaff10x). E. pustulosus and H. boettgeri chromosome-scale scaffolds were constructed with Dovetail Genomics Hi-C via the HiRise scaffolder173, followed by manual curation in Juicebox158,160,161 v1.9.0. Due to the fragmented nature of the E. coqui assembly, initial chromosome-scale scaffolds were rst constructed by synteny with E. pustulosus, then rened in Juicebox158,160,161 v1.9.0. Gaps in the E. pustulosus and H. boettgeri assemblies bridged by PacBio reads were resized using custom scripts (pbGapLen v0.0.2, https://bitbucket.org/rokhsar-lab/ xentr10/src/master/assembly) and lled with PBJelly174 (PBSuite v15.8.24). These two assemblies were polished with FreeBayes (v1.1.054-g49413aa) and ILEC (map4cns commit dd89f52, https://bitbucket. org/rokhsar-lab/map4cns). A nal round of gap-lling was then performed on the three assemblies using Platanus175 (v1.2.1). Previously published L. ailaonicum30 (GCA_018994145.1) and P. adspersus28 (GCA_004786255.1) assemblies were manually corrected in Juicebox158,160,161 (v1.11.08) using their respective Hi-C and Chicago data (Supplementary Data 1). Gaps in the corrected P. adspersus scaffolds were resized with PacBio reads (as described above) and lled using Platanus175 (v1.2.1) with published Illumina TruSeq PE data obtained from NCBI (PRJNA439445). As described elsewhere176, all assemblies were screened for contaminants before scaffolding, and only nal scaffolds and contigs longer than 1 kb were retained for downstream analyses. More details on assembly procedures can be found in (Supplementary Note 3). Genomic repeats in all ve species were annotated with RepeatMasker168,170 (v4.0.7 and v4.0.9) using the repeat library generated above. Protein-coding genes were annotated for E. coqui, E. pustulosus, H. boettgeri, and P. adspersus using the DOE-JGI IGC165 (v5.0) pipeline with homology and transcript evidence. For each respective species, newly generated RNA-seq data were combined with public H. boettgeri27 (BioProject PRJNA306175) and P. adspersus28 (BioProject PRJNA439445) data and E. coqui data (stages 7, 10, and 13 hindlimb [Harvard University]; stage 910 tail n skin [French National Center for Scientic Research]). TAs used as input to IGC were assembled with Trinity166,167 (v2.5.1) and ltered using the heuristics described in Supplementary Note 3. Synteny and ancestral chromosome inference One-to-one gene ortholog set between frog proteomes was obtained from the output from OrthoVenn264 (https://orthovenn2. bioinfotoolkits.net) using an E value of 1 105 and an ination value of 1.5 (Supplementary Note 4). The assemblies of all frog species and axolotl were pairwise aligned against the X. tropicalis genome using Cactus177 (commit e4d0859) (Supplementary Note 4). Pairwise collinear runs were merged into multiple sequence alignments with ROAST/MULTIZ178 (v012109) in order of phylogenetic topology from TimeTree179 (http://www.timetree.org), then sorted with LAST180 (v979) (Supplementary Note 4). 11 Article Phylogeny and estimation of sequence divergence Fourfold degenerate bases of one-to-one orthologs were obtained and reformatted from the MAFFT (v7.427) alignment as described in ref. 176 (Supplementary Note 4). The maximum-likelihood phylogeny was obtained with RAxML181 (v8.2.11) using the GTR+Gamma model of substitution with outgroup Ambystoma mexicanum. Divergence times were calculated with MEGA7182 (v7.0.26) with the GTR+Gamma model of substitution using Reltime method183. Chromosome evolution A custom script176 (cactus_lter.py v1.0, https://github.com/abmudd/ Assembly) was used to extract pairwise alignments from the ROASTmerged MAF le and convert alignments into runs of collinearity. The runs of collinearity were visualized with Circos184 (v0.69-6) (Supplementary Note 4) and JCVI185 (jcvi.graphics.karyotype v0.8.12, https:// github.com/tanghaibao/jcvi). Centromeres, satellites, and pericentromeric repeats Tandem repeats were called using Tandem Repeats Finder69 (v4.09; params: 2 5 7 80 10 50 2000 -l 6 -d -h -ngs). To identify tandem repeats enriched in pericentromeric and subtelomeric regions, we extracted the monomer sequences of all tandem repeats overlapping the region of interest. A database of non-redundant monomers was created by making a dimer database. Dimers were clustered with BlastClust186 v2.2.26 (-S 75 -p F -L 0.45 -b F -W 10). A non-redundant monomer database was created using the most common monomer size from each cluster. The non-redundant sequences were mapped to the genome with BLASTN187 (BLAST+ v2.9.0; -outfmt 6 -evalue 1e3). The enriched monomeric sequences in centromeres and subtelomeres were identied by selecting the highest normalized rations of tandem sequence footprints in the region of interest over the remaining portions of the genome. For more detail, see Supplementary Note 5. Genetic variation Reads were aligned with BWA-MEM188 (v0.7.17-r1188) and alignments were processed using SAMtools189 (v1.9-93-g0ca96a4), keeping only properly paired reads (samtools view -f3 -F3852) for variant calling. Variants were called with FreeBayes164 (v1.1.0-54-g49413aa; --standardlters --genotype-qualities --strict-vcf --report-monomorphic). Only biallelic SNPs with depth within mode 1.78SDs were retained. An allelebalance lter [0.30.7] for heterozygous genotypes was also applied. Segmental heterozygosity/homozygosity was estimated using windows of 500 kb with 50-kb step using BEDtools190 (v2.28.0) for pooled samples or snvrate191 (v2.0, https://bitbucket.org/rokhsar-lab/wgsanalysis). For more detail, see Supplementary Note 2. GC content, gene, and repeat landscape GC-content percentages were calculated in 1-Mb bins sliding every 50 kb. Gene densities were obtained using a window size of 250 kb sliding every 12.5 kb. The repeat density matrix for X. tropicalis was obtained by counting base pairs per 1 Mb (sliding every 200 kb) covered by repeat families and classes of repeats. The principal component analysis (PCA) was performed on the density matrix composed of 7253 overlapping 1-Mb bins and 3070 repeats (Supplementary Note 5). The rst (PC1) and second (PC2) components were smoothed using a cubic spline method. Chromatin immunoprecipitation Xenopus tropicalis XTN-6 cells192 were grown in 70% calcium-free L-15 (US Biologicals cat# L2101-02-50L), pH 7.2/10% Fetal Bovine Serum/ Penicillin-Streptomycin (Invitrogen cat# 15140-163) at RT. Native MNase ChIP-seq protocol was performed as described previously in Smith et al.88. Approximately 40 million cells were trypsinized and collected; nuclei were isolated by dounce extraction and collected with a sucrose cushion. Chromatin was digested to mononucleosomes by Nature Communications | (2024)15:579 https://doi.org/10.1038/s41467-023-43012-9 MNase. Nuclei were lysed and soluble nucleosomes were extracted overnight at 4 C. Extracted mononucleosomes were precleared with Protein A dynabeads (Invitrogen cat# 100-02D) for at least 4 h at 4 C. A sample was taken for input after pre-clearing. Protein A dynabeads were bound to 10-g antibody (50 g/L nal concentration of either Rb-anti-Xl Cenp-a [cross-reactive with X. tropicalis], Rb-anti-H4 Abcam cat# 7311, or Rb-anti-H3 Abcam cat# 1791) and incubated overnight with precleared soluble mononucleosomes at 4 C. Dynabeads bound to 50 g/L nal concentration of Rabbit IgG antibody (Jackson ImmunoResearch cat# 011-000-003) were collected with a magnet and washed three times with TBST (0.1% Triton X-100) before elution with 0.1% SDS in TE and proteinase K incubation at 65 C with shaking for at least 4 h. Isolated and input mononucleosomes were size-selected using Ampure beads (Beckman cat# A63880) and prepared for sequencing using the NEBNext Ultra II DNA Library Prep Kit for Illumina (NEB cat# E7654). Three replicates were sequenced on an Illumina HiSeq 4000 lane 2 150 bp by the Stanford Functional Genomics Facility. PE reads were trimmed with Trimmomatic193 (v0.39), removing universal Illumina primers and Nextera-PE indices. Processed PE reads were mapped with Minimap2194 (v2.17-r941) against the unmasked genome reference. SAMtools189 (v1.9-93-g0ca96a4) was used for sorting and indexing the alignments. Read counts (mapping quality [MQ] 0) per 10-kb bin (nonoverlapping) for all samples were calculated with multiBamSummary from deepTools195 (v3.3.0). Read counts were normalized by the total number of counts in the chromosomes per sample (Supplementary Note 5). Peaks were called with MACS2196 (v2.2.7.1) and custom scripts (https://bitbucket.org/rokhsarlab/xentr10/src/master/chipseq). Recombination and extended subtelomeres The reads from the F2 mapping population25 were aligned to the v10 genome sequence using BWA-MEM188 (v0.7.17-r1188). Variants were called using FreeBayes164 (v1.1.0-54-g49413aa; --standard-lters --genotype-qualities --strict-vcf ). SNPs were ltered, and valid F2 mapping sites were selected when the genotypes of the Nigerian F0 and the ICB F0 were xed and different and there was a depth of at least 10 for each F0 SNP. Maps were calculated using JoinMap197 v4.1 (Supplementary Note 5, Supplementary Data 2). The variation on the linkage map was smoothed using the not-a-knot cubic spline function calculated every 500 kb. The Pearson correlation coefcient, r, was calculated between recombination rates and genomic features that include GC content, repeat densities, and densities of reported CTCF and recombination hotspots198,199. Chromatin conformations and higher-order interactions Hi-C read pairs were mapped with Juicer158,159 (commit d3ee11b) and observed counts were extracted at 1 Mb resolution with Juicer Tools (commit d3ee11b). Centromeres were estimated manually in Juicebox160 and rened with Centurion200 v0.1.0-3-g985439c using ICE-balanced MQ 0 matrices (https://bitbucket.org/rokhsar-lab/xentr10/src/master/ hic). Rabl-like chromatin structure was visualized with PCA from KnightRuiz201-balanced MQ 30 matrices and signicance was estimated by permutation testing (10,000 iterations, one-sided = 0.01) using custom R202 scripts. Rabl-like constraint between p- and q-arms was measured as the sum of square distances (SSD) in PC1-PC2 dimensions, calculated between nonoverlapping bins traveling sequentially away from the centromere. Inter-/intra-chromosomal contact enrichment analyses were quantied from MQ 30 matrices using 2 tests in R v3.5.0 (hic-analysis.R v1.0, https://bitbucket.org/rokhsar-lab/xentr10/ src/master/hic). See Supplementary Note 5 for more details. A/B compartments A/B compartments were called with custom R202 scripts (call-compartments.R v0.1.0, https://bitbucket.org/bredeson/artisanal) from KnightRuiz-balanced (observed/expected normalized) MQ 30 Hi-C 12 Article contact correlation matrices generated with Juicer158,159 (Supplementary Note 5). Pearsons correlation between PC1 from the Hi-C correlation matrix and gene density was used to designate A and B compartments per chromosome. https://doi.org/10.1038/s41467-023-43012-9 8. 9. Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. 10. Data availability Data supporting the ndings of this work are available throughout the main text, Methods, Supplementary Information, Supplementary Data, or archived in Zenodo (https://doi.org/10.5281/zenodo. 8393403). All newly generated assemblies, annotations, and raw data are deposited in the NCBI GenBank and SRA databases: X. tropicalis under BioProject accession codes PRJNA577946 and PRJNA526297, E. coqui under BioProject accession code PRJNA578591, E. pustulosus under BioProject accession code PRJNA578590, and H. boettgeri under BioProject accession code PRJNA578589. L. ailaonicum and P. adspersus re-assemblies were deposited at NCBI GenBank under accession DAJOPU000000000 and DYDO00000000, respectively; the versions described in this manuscript are DAJOPU010000000 [https://www.ncbi.nlm.nih.gov/nuccore/ DAJOPU000000000.1] and DYDO01000000 [https://www.ncbi. nlm.nih.gov/nuccore/DYDO00000000.1]. Raw X. tropicalis ChIP-seq data are available at the NCBI SRA under BioProject accession code PRJNA726269 and the processed data via the NCBI GEO database under series accession GSE199671. The E. coqui tail n RNA-seq data generated in this study have been deposited in the NCBI SRA database under accession code PRJNA1022815. The E. coqui hindlimb developmental series RNA-seq data are available under restricted access as the project is not yet published, access can be obtained by contacting Mara Laslo at ml125@wellesley.edu. Source data are provided with this paper. 11. 12. 13. 14. 15. 16. 17. 18. 19. Code availability All custom scripts used in this work are archived203 in Zenodo at https:// doi.org/10.5281/zenodo.8393403 and can be found via the project repository at https://bitbucket.org/rokhsar-lab/xentr10 (tag v1.0) or via the individual repositories linked therein: https://github.com/abmudd/ Assembly, https://bitbucket.org/bredeson/artisanal, https://bitbucket. org/rokhsar-lab/map4cns, https://bitbucket.org/rokhsar-lab/wgsanalysis, https://bitbucket.org/rokhsar-lab/gbs-analysis, and https:// gitlab.com/Bredeson/wombat. 20. 21. 22. 23. References 1. 2. 3. 4. 5. 6. 7. Cannatella, D. C. & de S, R. O. Xenopus laevis as a model organism. Syst. Biol. 42, 476507 (1993). Beetschen, J. C. How did urodele embryos come into prominence as a model system? Int. J. Dev. Biol. 40, 629636 (1996). Brown, D. D. A tribute to the Xenopus laevis oocyte and egg. J. Biol. Chem. 279, 4529145299 (2004). Harland, R. M. & Grainger, R. M. Xenopus research: metamorphosed by genetics and genomics. Trends Genet. 27, 507515 (2011). Gurdon, J. B. & Hopwood, N. The introduction of Xenopus laevis into developmental biology: of empire, pregnancy testing and ribosomal genes. Int. J. Dev. Biol. 44, 4350 (2000). Blaustein, A. R. & Dobson, A. A message from the frogs. Nature 439, 143144 (2006). Farrer, R. A. et al. Multiple emergences of genetically diverse amphibian-infecting chytrids include a globalized hypervirulent recombinant lineage. Proc. Natl. Acad. Sci. USA 108, 1873218736 (2011). Nature Communications | (2024)15:579 24. 25. 26. 27. 28. 29. Whiles, M. R. et al. Disease-driven amphibian declines alter ecosystem processes in a tropical stream. Ecosystems 16, 146157 (2013). Gomes, A. et al. Bioactive molecules from amphibian skin: their biological activities with reference to therapeutic potentials for possible drug development. Indian J. Exp. Biol. 45, 579593 (2007). McCallum, M. L. Amphibian decline or extinction? Current declines dwarf background extinction rate. hpet 41, 483491 (2007). Ryan, M. J., Fox, J. H., Wilczynski, W. & Rand, A. S. Sexual selection for sensory exploitation in the frog Physalaemus pustulosus. Nature 343, 6667 (1990). Minsuk, S. B. & Keller, R. E. Surface mesoderm in Xenopus: a revision of the stage 10 fate map. Dev. Genes Evol. 207, 389401 (1997). Daczewska, M. & Saczko, J. Various DNA content in myotube nuclei during myotomal myogenesis in Hymenochirus boettgeri (Anura: Pipidae). Folia Biol. 51, 151157 (2003). Romero-Carvajal, A. et al. Embryogenesis and laboratory maintenance of the foam-nesting tngara frogs, genus Engystomops (= Physalaemus). Dev. Dyn. 238, 14441454 (2009). Ryan, M. J. The brain as a source of selection on the social niche: examples from the psychophysics of mate choice in tngara frogs. Integr. Comp. Biol. 51, 756770 (2011). Elinson, R. P. Metamorphosis in a frog that does not have a tadpole. Curr. Top. Dev. Biol. 103, 259276 (2013). Conlon, J. M. & Mechkarska, M. Host-defense peptides with therapeutic potential from skin secretions of frogs from the family Pipidae. Pharmaceuticals 7, 5877 (2014). Ryan, M. J. & Guerra, M. A. The mechanism of sound production in tngara frogs and its role in sexual selection and speciation. Curr. Opin. Neurobiol. 28, 5459 (2014). Womble, M., Pickett, M. & Nascone-Yoder, N. Frogs as integrative models for understanding digestive organ development and evolution. Semin. Cell Dev. Biol. 51, 92105 (2016). Burmeister, S. S. Neurobiology of female mate choice in frogs: auditory ltering and valuation. Integr. Comp. Biol. 57, 857864 (2017). Miller, K. E., Session, A. M. & Heald, R. Kif2a scales meiotic spindle size in Hymenochirus boettgeri. Curr. Biol. 29, 37203727.e5 (2019). Ferguson-Smith, M. A. & Trifonov, V. Mammalian karyotype evolution. Nat. Rev. Genet. 8, 950962 (2007). Zhang, G. et al. Comparative genomics reveals insights into avian genome evolution and adaptation. Science 346, 13111320 (2014). Kiazim, L. G. et al. Comparative mapping of the macrochromosomes of eight avian species provides further insight into their phylogenetic relationships and avian karyotype evolution. Cells 10, 362 (2021). Mitros, T. et al. A chromosome-scale genome assembly and dense genetic map for Xenopus tropicalis. Dev. Biol. 452, 820 (2019). Niu, L. et al. Three-dimensional folding dynamics of the Xenopus tropicalis genome. Nat. Genet. 53, 10751087 (2021). Session, A. M. et al. Genome evolution in the allotetraploid frog Xenopus laevis. Nature 538, 336343 (2016). Denton, R. D., Kudra, R. S., Malcom, J. W., Du Preez, L. & Malone, J. H. The African Bullfrog (Pyxicephalus adspersus) genome unites the two ancestral ingredients for making vertebrate sex chromosomes. Cold Spring Harb. Lab. 329847 https://doi.org/10.1101/ 329847 (2018). Li, J. et al. Genomic and transcriptomic insights into molecular basis of sexually dimorphic nuptial spines in Leptobrachium leishanense. Nat. Commun. 10, 5551 (2019). 13 Article 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. Li, Y. et al. Chromosome-level assembly of the mustache toad genome using third-generation DNA sequencing and Hi-C analysis. Gigascience 8, giz114 (2019). Lu, B. et al. A large genome with chromosome-scale assembly sheds light on the evolutionary success of a true toad (Bufo gargarizans). Mol. Ecol. Resour. 21, 12561273 (2021). Sun, Y.-B., Zhang, Y. & Wang, K. Perspectives on studying molecular adaptations of amphibians in the genomic era. Zool. Res 41, 351364 (2020). Wilson, A. C., Sarich, V. M. & Maxson, L. R. The importance of gene rearrangement in evolution: evidence from studies on rates of chromosomal, protein, and anatomical evolution. Proc. Natl. Acad. Sci. USA 71, 30283030 (1974). Gregory, T. R. Animal genome size database. http://www. genomesize.com (2023). Sotero-Caio, C. G., Challis, R., Kumar, S. & Blaxter, M. Genomes on a Tree (GoaT): a centralized resource for eukaryotic genome sequencing initiatives. BISS 5, e74138 (2021). Morescalchi, A. Evolution and karyology of the amphibians. Boll. Zool. 47, 113126 (1980). Bush, G. L., Case, S. M., Wilson, A. C. & Patton, J. L. Rapid speciation and chromosomal evolution in mammals. Proc. Natl. Acad. Sci. USA 74, 39423946 (1977). Nowoshilow, S. et al. The axolotl genome and the evolution of key tissue formation regulators. Nature 554, 5055 (2018). Smith, J. J. et al. A chromosome-scale assembly of the axolotl genome. Genome Res. 29, 317324 (2019). Nrnberger, B. et al. A dense linkage map for a large repetitive genome: discovery of the sex-determining region in hybridizing re-bellied toads (Bombina bombina and Bombina variegata). G3 11, jkab286 (2021). Deakin, J. E., Graves, J. A. M. & Rens, W. The evolution of marsupial and monotreme chromosomes. Cytogenet. Genome Res. 137, 113129 (2012). Damas, J. et al. Evolution of the ancestral mammalian karyotype and syntenic regions. Proc. Natl. Acad. Sci. USA 119, e2209139119 (2022). OConnor, R. E. et al. Reconstruction of the diapsid ancestral genome permits chromosome evolution tracing in avian and nonavian dinosaurs. Nat. Commun. 9, 1883 (2018). Bogart, J. P., Balon, E. K. & Bruton, M. N. The chromosomes of the living coelacanth and their remarkable similarity to those of one of the most ancient frogs. J. Hered. 85, 322325 (1994). Hellsten, U. et al. The genome of the Western clawed frog Xenopus tropicalis. Science 328, 633636 (2010). Carneiro, M. O. et al. Pacic biosciences sequencing technology for genotyping and variation discovery in human data. BMC Genomics 13, 375 (2012). Koren, S. et al. Hybrid error correction and de novo assembly of single-molecule sequencing reads. Nat. Biotechnol. 30, 693700 (2012). Quail, M. A. et al. A tale of three next generation sequencing platforms: comparison of Ion Torrent, Pacic Biosciences and Illumina MiSeq sequencers. BMC Genomics 13, 341 (2012). Loomis, E. W. et al. Sequencing the unsequenceable: expanded CGG-repeat alleles of the fragile X gene. Genome Res. 23, 121128 (2013). Feng, Y.-J. et al. Phylogenomics reveals rapid, simultaneous diversication of three major clades of Gondwanan frogs at the Cretaceous-Paleogene boundary. Proc. Natl. Acad. Sci. USA 114, E5864E5870 (2017). Schmid, M. et al. The chromosomes of Terraranan frogs. Insights into vertebrate cytogenetics. Cytogenet. Genome Res. 130, 114 (2010). Nature Communications | (2024)15:579 https://doi.org/10.1038/s41467-023-43012-9 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68. 69. 70. 71. 72. 73. 74. 75. 76. Rabello, M. N. Chromosomal studies in Brazilian anurans. Caryologia 23, 4559 (1970). Scheel, J. J. The chromosomes of some African anuran species. In Genetics and Mutagenesis of Fish (ed Schrder, J. H.) 113116 (Springer, Berlin, Heidelberg, 1973). Mezzasalma, M., Glaw, F., Odierna, G., Petraccioli, A. & Guarino, F. M. Karyological analyses of Pseudhymenochirus merlini and Hymenochirus boettgeri provide new insights into the chromosome evolution in the anuran family Pipidae. Zoologischer Anz.A J. Comp. Zool. 258, 4753 (2015). Temple, G. et al. The completion of the mammalian gene collection (MGC). Genome Res. 19, 23242333 (2009). Marin, R. et al. Convergent origination of a Drosophila-like dosage compensation mechanism in a reptile lineage. Genome Res. 27, 19741987 (2017). Owens, N. D. L. et al. Measuring absolute RNA copy numbers at high temporal resolution reveals transcriptome kinetics in development. Cell Rep. 14, 632647 (2016). Warren, W. C. et al. A new chicken genome assembly provides insight into avian genome structure. G3 7, 109117 (2017). Howe, K. et al. The zebrash reference genome sequence and its relationship to the human genome. Nature 496, 498503 (2013). Mouse Genome Sequencing Consortium. Initial sequencing and comparative analysis of the mouse genome. Nature 420, 520562 (2002). Lander, E. S. et al. Initial sequencing and analysis of the human genome. Nature 409, 860921 (2001). Venter, J. C. et al. The sequence of the human genome. Science 291, 13041351 (2001). Lovell, P. V. et al. Conserved syntenic clusters of protein coding genes are missing in birds. Genome Biol. 15, 565 (2014). Xu, L. et al. OrthoVenn2: a web server for whole-genome comparison and annotation of orthologous clusters across multiple species. Nucleic Acids Res. 47, W52W58 (2019). Hartley, G. & ONeill, R. Centromere repeats: Hidden gems of the genome. Genes 10, 223 (2019). Chueh, A. C., Wong, L. H., Wong, N. & Choo, K. H. A. Variable and hierarchical size distribution of L1-retroelement-enriched CENP-A clusters within a functional human neocentromere. Hum. Mol. Genet. 14, 8593 (2005). Kuznetsova, I. S. et al. LINE-related component of mouse heterochromatin and complex chromocenters composition. Chromosome Res. 24, 309323 (2016). Suh, A. The specic requirements for CR1 retrotransposition explain the scarcity of retrogenes in birds. J. Mol. Evol. 81, 1820 (2015). Benson, G. Tandem Repeats Finder: A program to analyze DNA sequences. Nucleic Acids Res. 27, 573580 (1999). Nagylaki, T. Introduction to Theoretical Population Genetics (Springer Berlin Heidelberg, 1992). Igawa, T. et al. Inbreeding ratio and genetic relationships among strains of the Western clawed frog, Xenopus tropicalis. PLoS ONE 10, e0133963 (2015). Ford, L. S. & Cannatella, D. C. The major clades of frogs. Herpetol. Monogr. 7, 94117 (1993). Bhutkar, A. et al. Chromosomal rearrangement inferred from comparisons of 12 Drosophila genomes. Genetics 179, 16571680 (2008). Pyron, R. A. Divergence time estimation using fossils as terminal taxa and the origins of Lissamphibia. Syst. Biol. 60, 466481 (2011). Wright, S. On the probability of xation of reciprocal translocations. Am. Nat. 75, 513522 (1941). Lande, R. The xation of chromosomal rearrangements in a subdivided population with local extinction and colonization. Heredity 54, 323332 (1985). 14 Article 77. 78. 79. 80. 81. 82. 83. 84. 85. 86. 87. 88. 89. 90. 91. 92. 93. 94. 95. 96. 97. 98. Schubert, I. & Lysak, M. A. Interpretation of karyotype evolution should consider chromosome structural constraints. Trends Genet. 27, 207216 (2011). Lysak, M. A. Celebrating Mendel, McClintock, and Darlington: on end-to-end chromosome fusions and nested chromosome fusions. Plant Cell 34, 24752491 (2022). Grifn, D. K., Robertson, L. B. W., Tempest, H. G. & Skinner, B. M. The evolution of the avian genome as revealed by comparative molecular cytogenetics. Cytogenet. Genome Res. 117, 6477 (2007). Deakin, J. E. & Ezaz, T. Understanding the evolution of reptile chromosomes through applications of combined cytogenetics and genomics approaches. Cytogenet. Genome Res. 157, 720 (2019). Maruyama, T. & Imai, H. T. Evolutionary rate of the mammalian karyotype. J. Theor. Biol. 90, 111121 (1981). Olmo, E. Rate of chromosome changes and speciation in reptiles. Genetica 125, 185203 (2005). Duret, L. & Galtier, N. Biased gene conversion and the evolution of mammalian genomic landscapes. Annu. Rev. Genomics Hum. Genet. 10, 285311 (2009). Bogart, J. P. The Inuence of Life History on Karyotypic Evolution in Frogs (Academic Press, Inc., 1991). Bogart, J. P. & Hedges, S. B. Rapid chromosome evolution in Jamaican frogs of the genus Eleutherodactylus (Leptodactylidae). J. Zool. 235, 931 (1995). Jagannathan, M., Cummings, R. & Yamashita, Y. M. A conserved function for pericentromeric satellite DNA. eLife 7, e34122 (2018). Edwards, N. S. & Murray, A. W. Identication of Xenopus CENP-A and an associated centromeric DNA repeat. Mol. Biol. Cell 16, 18001810 (2005). Smith, O. K. et al. Identication and characterization of centromeric sequences in Xenopus laevis. Cold Spring Harb. Lab. https://doi.org/10.1101/2020.06.23.167643 (2020). Penke, T. J. R., McKay, D. J., Strahl, B. D., Matera, A. G. & Duronio, R. J. Direct interrogation of the role of H3K9 in metazoan heterochromatin function. Genes Dev. 30, 18661880 (2016). Di Giacomo, M. et al. Multiple epigenetic mechanisms and the piRNA pathway enforce LINE1 silencing during adult spermatogenesis. Mol. Cell 50, 601608 (2013). Drau, A., Venu, V., Avdievich, E., Gaspar, L. & Jones, F. C. Genome-wide recombination map construction from single individuals using linked-read sequencing. Nat. Commun. 10, 4309 (2019). Backstrom, N. et al. The recombination landscape of the zebra nch Taeniopygia guttata genome. Genome Res. 20, 485495 (2010). Groenen, M. A. M. et al. A high-density SNP-based linkage map of the chicken genome reveals sequence features correlated with recombination rate. Genome Res. 19, 510519 (2009). Schield, D. R. et al. Snake recombination landscapes are concentrated in functional regions despite PRDM9. Mol. Biol. Evol. 37, 12721294 (2020). Kong, A. et al. A high-resolution recombination map of the human genome. Nat. Genet. 31, 241247 (2002). Campbell, C. L., Bhrer, C., Morrow, B. E., Boyko, A. R. & Auton, A. A pedigree-based map of recombination in the domestic dog genome. G3 6, 35173524 (2016). Tortereau, F. et al. A high density recombination map of the pig reveals a correlation between sex-specic recombination and GC content. BMC Genomics 13, 586 (2012). Jensen-Seaman, M. I. et al. Comparative recombination rates in the rat, mouse, and human genomes. Genome Res. 14, 528538 (2004). Nature Communications | (2024)15:579 https://doi.org/10.1038/s41467-023-43012-9 99. 100. 101. 102. 103. 104. 105. 106. 107. 108. 109. 110. 111. 112. 113. 114. 115. 116. 117. 118. 119. 120. Baker, Z. et al. Repeated losses of PRDM9-directed recombination despite the conservation of PRDM9 across vertebrates. eLife 6, e24133 (2017). Kuhl, L.-M. & Vader, G. Kinetochores, cohesin, and DNA breaks: Controlling meiotic recombination within pericentromeres. Yeast 36, 121127 (2019). Termolino, P., Cremona, G., Consiglio, M. F. & Conicella, C. Insights into epigenetic landscape of recombination-free regions. Chromosoma 125, 301308 (2016). Singhal, S. et al. Stable recombination hotspots in birds. Science 350, 928932 (2015). Galtier, N., Piganeau, G., Mouchiroud, D. & Duret, L. GC-content evolution in mammalian genomes: the biased gene conversion hypothesis. Genetics 159, 907911 (2001). Meunier, J. & Duret, L. Recombination drives the evolution of GCcontent in the human genome. Mol. Biol. Evol. 21, 984990 (2004). Lam, B. S. & Carroll, D. Tandemly repeated DNA sequences from Xenopus laevis. I. Studies on sequence organization and variation in satellite 1 DNA (741 base-pair repeat). J. Mol. Biol. 165, 567585 (1983). Cohen, S., Menut, S. & Mchali, M. Regulated formation of extrachromosomal circular DNA molecules during development in Xenopus laevis. Mol. Cell. Biol. 19, 66826689 (1999). Ogiwara, I. V.-S. I. N. Es A new superfamily of vertebrate SINEs that are widespread in vertebrate genomes and retain a strongly conserved segment within each repetitive unit. Genome Res. 12, 316324 (2002). Rao, S. S. P. et al. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell 159, 16651680 (2014). Mascher, M. et al. A chromosome conformation capture ordered sequence of the barley genome. Nature 544, 427433 (2017). Hoencamp, C. et al. 3D genomics across the tree of life reveals condensin II as a determinant of architecture type. Science 372, 984989 (2021). Rabl, C. ber Zelltheilung. Morphologisches Jahrbuch 10, 214330 (1885). Muller, H., Gil, J. Jr & Drinnenberg, I. A. The impact of centromeres on spatial genome architecture. Trends Genet. 35, 565578 (2019). Sperling, K. & Ldtke, E. K. Arrangement of prematurely condensed chromosomes in cultured cells and lymphocytes of the Indian muntjac. Chromosoma 83, 541553 (1981). Cremer, T. et al. Rabls model of the interphase chromosome arrangement tested in Chinese hamster cells by premature chromosome condensation and laser-UV-microbeam experiments. Hum. Genet. 60, 4656 (1982). Mathog, D., Hochstrasser, M., Gruenbaum, Y., Saumweber, H. & Sedat, J. Characteristic folding pattern of polytene chromosomes in Drosophila salivary gland nuclei. Nature 308, 414421 (1984). Stevens, T. J. et al. 3D structures of individual mammalian genomes studied by single-cell Hi-C. Nature 544, 5964 (2017). Dudchenko, O. et al. De novo assembly of the Aedes aegypti genome using Hi-C yields chromosome-length scaffolds. Science 356, 9295 (2017). Funabiki, H., Hagan, I., Uzawa, S. & Yanagida, M. Cell cycledependent specic positioning and clustering of centromeres and telomeres in ssion yeast. J. Cell Biol. 121, 961976 (1993). Duan, Z. et al. A three-dimensional model of the yeast genome. Nature 465, 363367 (2010). Armstrong, S. J., Franklin, F. C. & Jones, G. H. Nucleolusassociated telomere clustering and pairing precede meiotic chromosome synapsis in Arabidopsis thaliana. J. Cell Sci. 114, 42074217 (2001). 15 Article 121. 122. 123. 124. 125. 126. 127. 128. 129. 130. 131. 132. 133. 134. 135. 136. 137. 138. 139. 140. 141. Santos, A. P. & Shaw, P. Interphase chromosomes and the Rabl conguration: does genome size matter? J. Microsc. 214, 201206 (2004). Cowan, C. R., Carlton, P. M. & Cande, W. Z. The polar arrangement of telomeres in interphase and meiosis. Rabl organization and the bouquet. Plant Physiol. 125, 532538 (2001). Therizols, P., Duong, T., Dujon, B., Zimmer, C. & Fabre, E. Chromosome arm length and nuclear constraints determine the dynamic relationship of yeast subtelomeres. Proc. Natl. Acad. Sci. USA 107, 20252030 (2010). Buttrick, G. J. et al. Nsk1 ensures accurate chromosome segregation by promoting association of kinetochores to spindle poles during anaphase B. Mol. Biol. Cell 22, 44864502 (2011). Dernburg, A. F. et al. Perturbation of nuclear architecture by longdistance chromosome interactions. Cell 85, 745759 (1996). Hiraoka, Y. et al. The onset of homologous chromosome pairing during Drosophila melanogaster embryogenesis. J. Cell Biol. 120, 591600 (1993). Marshall, W. F., Dernburg, A. F., Harmon, B., Agard, D. A. & Sedat, J. W. Specic interactions of chromatin with the nuclear envelope: positional determination within the nucleus in Drosophila melanogaster. Mol. Biol. Cell 7, 825842 (1996). Rowley, M. J. & Corces, V. G. Organizational principles of 3D genome architecture. Nat. Rev. Genet. 19, 789800 (2018). Lu, J. Y. et al. Homotypic clustering of L1 and B1/Alu repeats compartmentalizes the 3D genome. Cell Res. 31, 613630 (2021). Fishman, V. et al. 3D organization of chicken genome demonstrates evolutionary conservation of topologically associated domains and highlights unique architecture of erythrocytes chromatin. Nucleic Acids Res. 47, 648665 (2019). Kaaij, L. J. T., van der Weide, R. H., Ketting, R. F. & de Wit, E. Systemic loss and gain of chromatin architecture throughout zebrash development. Cell Rep. 24, 110.e4 (2018). Eagen, K. P., Aiden, E. L. & Kornberg, R. D. Polycomb-mediated chromatin loops revealed by a subkilobase-resolution chromatin interaction map. Proc. Natl. Acad. Sci. USA 114, 87648769 (2017). Dong, P. et al. 3D chromatin architecture of large plant genomes determined by local A/B compartments. Mol. Plant 10, 14971509 (2017). Francke, U. 2012 William Allan Award: adventures in cytogenetics. Am. J. Hum. Genet. 92, 325337 (2013). Uno, Y. et al. Diversity in the origins of sex chromosomes in anurans inferred from comparative mapping of sexual differentiation genes for three species of the Raninae and Xenopodinae. Chromosome Res. 16, 9991011 (2008). Uno, Y. et al. Inference of the protokaryotypes of amniotes and tetrapods and the evolutionary processes of microchromosomes from comparative gene mapping. PLoS ONE 7, e53027 (2012). Parada, L. A., McQueen, P. G., Munson, P. J. & Misteli, T. Conservation of relative chromosome positioning in normal and cancer cells. Curr. Biol. 12, 16921697 (2002). Parada, L. A., McQueen, P. G. & Misteli, T. Tissue-specic spatial organization of genomes. Genome Biol. 5, R44 (2004). Lieberman-Aiden, E. et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science 326, 289293 (2009). Uno, Y., Nishida, C., Takagi, C., Ueno, N. & Matsuda, Y. Homoeologous chromosomes of Xenopus laevis are highly conserved after whole-genome duplication. Heredity 111, 430436 (2013). Kozubek, S. et al. The topological organization of chromosomes 9 and 22 in cell nuclei has a determinative role in the induction of t(9,22) translocations and in the pathogenesis of t(9,22) leukemias. Chromosoma 108, 426435 (1999). Nature Communications | (2024)15:579 https://doi.org/10.1038/s41467-023-43012-9 142. Branco, M. R. & Pombo, A. Intermingling of chromosome territories in interphase suggests role in translocations and transcription-dependent associations. PLoS Biol. 4, e138 (2006). 143. Rosin, L. F. et al. Chromosome territory formation attenuates the translocation potential of cells. eLife 8, e49553 (2019). 144. Bright, A. R. et al. Combinatorial transcription factor activities on open chromatin induce embryonic heterogeneity in vertebrates. EMBO J. 40, e104913 (2021). 145. Kakebeen, A. D., Chitsazan, A. D., Williams, M. C., Saunders, L. M. & Wills, A. E. Chromatin accessibility dynamics and single cell RNASeq reveal new regulators of regeneration in neural progenitors. eLife 9, e52648 (2020). 146. del Pino, E. M. et al. A comparative analysis of frog early development. Proc. Natl. Acad. Sci. USA 104, 1188211888 (2007). 147. Vargas, A. & Del Pino, E. M. Analysis of cell size in the gastrula of ten frog species reveals a correlation of egg with cell sizes, and a conserved pattern of small cells in the marginal zone. J. Exp. Zool. B Mol. Dev. Evol. 328, 8896 (2017). 148. Oswald, P. et al. Locality, time and heterozygosity affect chytrid infection in yellow-bellied toads. Dis. Aquat. Organ. 142, 225237 (2020). 149. Alford, R. A., Dixon, P. M. & Pechmann, J. H. Ecology. Global amphibian population declines. Nature 412, 499500 (2001). 150. Leung, B. et al. Clustered versus catastrophic global vertebrate declines. Nature 588, 267271 (2020). 151. Gvodk, V., Knytl, M., Zassi-Boulou, A-G, Fornaini, N. R. & Bergelov, B. Tetraploidy in the Boettgers dwarf clawed frog (Pipidae: Hymenochirus boettgeri) from the Congo indicates non-conspecicity with the captive population, Zoological Journal of the Linnean Society zlad119 https://doi.org/10.1093/zoolinnean/ zlad119 (2023). 152. Weisenfeld, N. I., Kumar, V., Shah, P., Church, D. M. & Jaffe, D. B. Direct determination of diploid genome sequences. Genome Res. 27, 757767 (2017). 153. Ye, C., Hill, C. M., Wu, S., Ruan, J. & Ma, Z. S. DBG2OLC: Efcient assembly of large genomes using long erroneous reads of the third generation sequencing technologies. Sci. Rep. 6, 31900 (2016). 154. Koren, S. et al. Canu: Scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res. 27, 722736 (2017). 155. Kurtz, S. et al. Versatile and open software for comparing large genomes. Genome Biol. 5, R12 (2004). 156. Chakraborty, M., Baldwin-Brown, J. G., Long, A. D. & Emerson, J. J. Contiguous and accurate de novo assembly of metazoan genomes with modest long read coverage. Nucleic Acids Res. 44, e147e147 (2016). 157. Boetzer, M., Henkel, C. V., Jansen, H. J., Butler, D. & Pirovano, W. Scaffolding pre-assembled contigs using SSPACE. Bioinformatics 27, 578579 (2011). 158. Durand, N. C. et al. Juicer provides a one-click system for analyzing loop-resolution Hi-C experiments. Cell Syst. 3, 9598 (2016). 159. Tange, O. GNU Parallel 2018 (Lulu.com, 2018). 160. Durand, N. C. et al. Juicebox provides a visualization system for HiC contact maps with unlimited zoom. Cell Syst. 3, 99101 (2016). 161. Dudchenko, O., Shamim, M. S., Batra, S. S. & Durand, N. C. The Juicebox Assembly Tools module facilitates de novo assembly of mammalian genomes with chromosome-length scaffolds for under $1000. Preprint at https://www.biorxiv.org/content/10. 1101/254797v1 (2018). 162. Chin, C.-S. et al. Nonhybrid, nished microbial genome assemblies from long-read SMRT sequencing data. Nat. Methods 10, 563569 (2013). 16 Article 163. Walker, B. J. et al. Pilon: An integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS ONE 9, e112963 (2014). 164. Garrison, E. & Marth, G. Haplotype-based variant detection from short-read sequencing. Preprint at https://arxiv.org/abs/1207. 3907 (2012). 165. Shu, S., Rokhsar, D., Goodstein, D., Hayes, D. & Mitros, T. JGI Plant Genomics Gene Annotation Pipeline. https://www.osti.gov/biblio/ 1241222 (2014). 166. Grabherr, M. G. et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644652 (2011). 167. Haas, B. J. et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat. Protoc. 8, 14941512 (2013). 168. Smit, A. F. A. & Hubley, R. RepeatModeler Open-1.0. https://www. repeatmasker.org/RepeatModeler (20082015). 169. Jurka, J. et al. Repbase update, a database of eukaryotic repetitive elements. Cytogenet. Genome Res. 110, 462467 (2005). 170. Smit, A. F. A., Hubley, R. & Green, P. RepeatMasker Open-4.0. http://www.repeatmasker.org (20132015). 171. Chapman, J. A. et al. Meraculous: de novo genome assembly with short paired-end reads. PLoS ONE 6, e23501 (2011). 172. Goltsman, E., Ho, I. & Rokhsar, D. Meraculous-2D: haplotypesensitive assembly of highly heterozygous genomes. Preprint at https://arxiv.org/ftp/arxiv/papers/1703/1703.09852.pdf (2017). 173. Putnam, N. H. et al. Chromosome-scale shotgun assembly using an in vitro method for long-range linkage. Genome Res. 26, 342350 (2016). 174. English, A. C. et al. Mind the gap: upgrading genomes with Pacic Biosciences RS long-read sequencing technology. PLoS ONE 7, e47768 (2012). 175. Kajitani, R. et al. Efcient de novo assembly of highly heterozygous genomes from whole-genome shotgun short reads. Genome Res. 24, 13841395 (2014). 176. Mudd, A. B., Bredeson, J. V., Baum, R., Hockemeyer, D. & Rokhsar, D. S. Analysis of muntjac deer genome and chromatin architecture reveals rapid karyotype evolution. Commun. Biol. 3, 110 (2020). 177. Paten, B. et al. Cactus: Algorithms for genome multiple sequence alignment. Genome Res. 21, 15121528 (2011). 178. Blanchette, M. et al. Aligning multiple genomic sequences with the threaded blockset aligner. Genome Res. 14, 708715 (2004). 179. Kumar, S., Stecher, G., Suleski, M. & Hedges, S. B. TimeTree: a resource for timelines, timetrees, and divergence times. Mol. Biol. Evol. 34, 18121819 (2017). 180. Kiebasa, S. M., Wan, R., Sato, K., Horton, P. & Frith, M. C. Adaptive seeds tame genomic sequence comparison. Genome Res. 21, 487493 (2011). 181. Stamatakis, A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 13121313 (2014). 182. Kumar, S., Stecher, G. & Tamura, K. MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 18701874 (2016). 183. Tamura, K. et al. Estimating divergence times in large molecular phylogenies. Proc. Natl. Acad. Sci. USA 109, 1933319338 (2012). 184. Krzywinski, M. et al. Circos: an information aesthetic for comparative genomics. Genome Res. 19, 16391645 (2009). 185. Tang, H. et al. Synteny and collinearity in plant genomes. Science 320, 486488 (2008). 186. Dondoshansky, I. & Wolf, Y. Blastclust (NCBI Software Development Toolkit). ScienceOpen https://www.scienceopen.com/ document?vid=b654ab9a-231d-410a-832d-37c7c7bc7165 (2002). 187. Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinforma. 10, 421 (2009). Nature Communications | (2024)15:579 https://doi.org/10.1038/s41467-023-43012-9 188. Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. Preprint at https://arxiv.org/abs/1303. 3997 (2013). 189. Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 20782079 (2009). 190. Quinlan, A. R. BEDTools: the Swiss-army tool for genome feature analysis. Curr. Protoc. Bioinforma. 47, 11.12.134 (2014). 191. Bredeson, J. V. et al. Sequencing wild and cultivated cassava and related species reveals extensive interspecic hybridization and genetic diversity. Nat. Biotechnol. 34, 562570 (2016). 192. Gorbsky, G. J. et al. Developing immortal cell lines from Xenopus embryos, four novel cell lines derived from Xenopus tropicalis. Open Biol. 12, 19 (2022). 193. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a exible trimmer for Illumina sequence data. Bioinformatics 30, 21142120 (2014). 194. Li, H. Minimap2: Pairwise alignment for nucleotide sequences. Bioinformatics 34, 30943100 (2018). 195. Ramrez, F. et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 44, W160W165 (2016). 196. Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, R137 (2008). 197. Van Ooijen, J. W. Multipoint maximum likelihood mapping in a fullsib family of an outbreeding species. Genet. Res. 93, 343349 (2011). 198. Myers, S., Bottolo, L., Freeman, C., McVean, G. & Donnelly, P. A ne-scale map of recombination rates and hotspots across the human genome. Science 310, 321324 (2005). 199. Shifman, S. et al. A high-resolution single nucleotide polymorphism genetic map of the mouse genome. PLoS Biol. 4, e395 (2006). 200. Varoquaux, N. et al. Accurate identication of centromere locations in yeast genomes using Hi-C. Nucleic Acids Res. 43, 53315339 (2015). 201. Knight, P. A. & Ruiz, D. A fast algorithm for matrix balancing. IMA J. Numer. Anal. 33, 10291047 (2012). 202. R Core Team. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna Austria http://www.R-project.org/ (2013). 203. Bredeson, J. V. et al. Conserved chromatin and repetitive patterns reveal slow genome evolution in frogs. https://doi.org/10.5281/ zenodo.8393403 (2023). Acknowledgements We thank Karen Lundy and the Functional Genomics Laboratory at the University of California Berkeley for running quality control on extracted DNA and RNA and for preparing Illumina short-insert libraries; Oanh Nguyen and the DNA Technologies and Expression Analysis Cores at the University of California Davis Genome Center for preparing and sequencing PacBio libraries; Dovetail Genomics for providing the Hi-C library preparation kit, running quality control on Hi-C libraries, and preparing and sequencing Hi-C libraries; Shana McDevitt and the Vincent J. Coates Genomics Sequencing Laboratory at the University of California Berkeley for sequencing Hi-C and Illumina short-insert libraries; Shengqiang Shu for advice on the use of the IGC annotation pipeline. We thank Rick Elinson for providing E. coqui frogs and tissues. We thank Gary Gorbsky from the Oklahoma Medical Research Foundation and Marko Horb and the National Xenopus Resource at the MBL for providing the XTN-6 cell lines. We also thank Chunhui Hou and colleagues for permission to access their Hi-C data before publication. This study was supported by NIH grants R01HD080708 to D.S.R.; R01GM086321, R01HD065705 to D.S.R. and R.M.H.; R35GM127069 to R.M.H.; R35 GM118183 to R.H. A.B.M. was supported by NIH grants T32GM007127 and T32HG000047 and a David L. Boren Fellowship. D.S.R. is grateful for support from the Marthella Foskett Brown Chair in Biological Sciences; R.M.H., the C.H. Li 17 Article Distinguished Chair in Molecular and Cell Biology; and R.H., the Flora Lamson Hewlett chair in biochemistry. A.F.S. and O.K.S. were supported by R01GM074728, O.K.S. by NIH T32 GM113854-02 and NSF GRFP; M.K.K. and M.Lane by R01HD102186; J.H. by NSF grants DEB-1701591 and DBI1702263; M.Laslo, a Graduate Women in Science Fellowship; T.K. by the Basic Science Research Program, National Research Foundation of Korea (NRF), Ministry of Education (2018R1A6A1A03025810), Future-leading Project Research Fund (1.200094.01) of UNIST and the Institute for Basic Science (IBS-R022-D1); J.B.W. and H.S.P. by R01GM104853, R01HD085901; M.J.R. by NSF IOS-0910112; Smithsonian Tropical Research Institute; Clark Hubbs Regents Professorship; L.M.S. by the Centre National de la Recherche Scientique (PEPS ExoMod Triton) m National dHistoire Naturelle (Action Transversale du and the Museu Museum Cycles biologiques: Evolution et adaptation) and a Scientic council post-doctoral position to G.K. This work used the Vincent J. Coates Genomics Sequencing Laboratory at the University of California Berkeley, supported by NIH grant S10OD018174, and the DNA Technologies and Expression Analysis Cores at the University of California Davis Genome Center, supported by NIH grant S10OD010786. This research used the National Energy Research Scientic Computing Center, a Department of Energy Ofce of Science User Facility supported by contract number DE-AC02-05CH11231. L.M.S. acknowledges the Ecole Normale Superieure de PARIS genomic platform for RNA sequencing m and the PCIA high-performance computing platform at Museu National dHistoire Naturelle. https://doi.org/10.1038/s41467-023-43012-9 sequenced 10x Genomics, PacBio, and Illumina mate-pair libraries. D.H. prepared Hi-C libraries. R.D.D. and J.H.M. provided early access to the Pad assembly. N.B. (Eco) provided bioinformatic support. L.M.S. led the Eco efforts. R.M.H. and D.S.R. led the project. Competing interests D.S.R. is a member of the Scientic Advisory Board of, and a minor shareholder in, Dovetail Genomics LLC, which provides as a service the high-throughput chromatin conformation capture (Hi-C) technology used in this study. M.K.K. is President and co-founder of Victory Genomics, Inc. The remaining authors declare no competing interests. Additional information Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-43012-9. Correspondence and requests for materials should be addressed to Daniel S. Rokhsar. Peer review information Nature Communications thanks Mark Blaxter and Amy Sater for their contribution to the peer review of this work. A peer review le is available. Reprints and permissions information is available at http://www.nature.com/reprints Author contributions J.V.B., A.B.M., S.M.R., T.M., R.M.H. and D.S.R. wrote the manuscript with feedback from M.Laslo, H.P.S., J.H., J.B.L., J.B.W., M.J.R., O.K.S., D.R.B., M.G.P., J.H., N.B., T.K., L.M.S., R.H., J.S., M.K.K., A.F.S. and D.H. Genomes were assembled by J.V.B., S.S.B. (Xtr); A.B.M., and K.C.B. (other frogs). S.M.R., A.B.M. and G.K. assembled transcripts and annotated genomes. S.M.R. and J.V.B. assessed gene completeness; S.M.R. analyzed repeat and recombination landscapes. S.M.R. and J.P. identied centromeric repeats. O.K.S., G.A.F. and A.F.S. conducted ChIP-seq experiments, and S.M.R. performed analysis. J.V.B. analyzed Hi-C features. T.M. constructed the linkage map. T.M. and J.V.B. analyzed heterozygosity. A.B.M. performed genome-wide comparisons. K.E.M. and R.H. examined Hbo metaphase spreads. M.K.K. and M.Lane inbred Xtr frogs. R.M.H. (Xtr); M.G.P. (Epu); K.E.M. and R.H. (Hbo); M.Laslo and J.H. (Eco) collected frogs. R.M.H. (Xtr); M.G.P., H.S.P. (Epu); and D.R.B. (Eco) collected tissue samples. A.B.M., D.R.B. (Eco); J.B.L. and I.P. (Xtr) extracted DNA. A.B.M., S.M.R. (Epu); K.E.M., R.H. (Hbo); and L.M.S. (Eco) extracted RNA and libraries were prepared by A.B.M. (Epu). M.Laslo, J.H. (Eco); K.E.M. and R.H. (Hbo) provided RNA-seq data. T.K., M.J.R., J.B.W. (Epu); and J.B.L. (Xtr) coordinated sequencing. C.P., J.G. and J.S. prepared and Publishers note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional afliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the articles Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the articles Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/. The Author(s) 2024 1 Department of Molecular and Cell Biology, Weill Hall, University of California, Berkeley, CA 94720, USA. 2DOE-Joint Genome Institute, 1 Cyclotron Road, Berkeley, CA 94720, USA. 3Department of Biochemistry, Stanford University School of Medicine, 279 Campus Drive, Beckman Center 409, Stanford, CA 94305-5307, USA. 4Computer Science Division, University of California Berkeley, 2626 Hearst Avenue, Berkeley, CA 94720, USA. 5HudsonAlpha Genome Sequencing Center, HudsonAlpha Institute for Biotechnology, Huntsville, AL 35806, USA. 6Pediatric Genomics Discovery Program, Departments of Pediatrics and Genetics, Yale University School of Medicine, 333 Cedar Street, New Haven, CT 06510, USA. 7Department of Organismic and Evolutionary Biology, and Museum of Comparative Zoology, Harvard University, Cambridge, MA 02138, USA. 8Dpartement Adaptation du Vivant, UMR 7221 CNRS, Musum National dHistoire Naturelle, Paris, France. 9Department of Biological Sciences, University of Cincinnati, Cincinnati, OH, USA. 10Department of Biomedical Engineering, Ulsan National Institute of Science and Technology, Ulsan 44919, Republic of Korea. 11Center for Genomic Integrity, Institute for Basic Science (IBS), Ulsan 44919, Republic of Korea. 12Department of Integrative Biology, Patterson Labs, 2401 Speedway, University of Texas, Austin, TX 78712, USA. 13Department of Biological Sciences, University of the Pacic, 3601 Pacic Avenue, Stockton, CA 95211, USA. 14Department of Molecular and Cell Biology and Institute of Systems Genomics, University of Connecticut, 181 Auditorium Road, Unit 3197, Storrs, CT 06269, USA. 15Department of Molecular Biosciences, Patterson Labs, 2401 Speedway, The University of Texas at Austin, Austin, TX 78712, USA. 16Innovative Genomics Institute, University of California, Berkeley, CA 94720, USA. 17Chan-Zuckerberg BioHub, 499 Illinois Street, San Francisco, CA 94158, USA. 18Okinawa Institute of Science and Technology Graduate University, Onna, Okinawa 9040495, Japan. 19These authors contributed equally: Jessen V. Bredeson, Austin B. Mudd, Soa Medina-Ruiz. e-mail: dsrokhsar@gmail.com Nature Communications | (2024)15:579 18  ...

... The item referenced in this repository content can be found by following the link on the descriptive page.  ...

... Breast Cancer Research Edwards et al. Breast Cancer Research (2024) 26:34 https://doi.org/10.1186/s13058-024-01791-z Open Access RESEARCH PTHrP intracrine actions divergently influence breast cancer growth through p27 and LIFR Courtney M. Edwards1,2, Jeremy F. Kane1,2, Jailyn A. Smith2,3, Dja M. Grant1,2,4, Jasmine A. Johnson2,3, Maria A. Hernandez Diaz2,3, Lawrence A. Vecchi III2,3, Kai M. Bracey5, Tolu N. Omokehinde1,2, Joseph R. Fontana2,6, Breelyn A. Karno2,6, Halee T. Scott2,6, Carolina J. Vogel7,8, Jonathan W. Lowery7,8,9,10, T. John Martin11,12 and Rachelle W. Johnson1,2,3* Abstract The role of parathyroid hormone (PTH)-related protein (PTHrP) in breast cancer remains controversial, with reports of PTHrP inhibiting or promoting primary tumor growth in preclinical studies. Here, we provide insight into these conflicting findings by assessing the role of specific biological domains of PTHrP in tumor progression through stable expression of PTHrP (-36-139aa) or truncated forms with deletion of the nuclear localization sequence (NLS) alone or in combination with the C-terminus. Although the full-length PTHrP molecule (-36-139aa) did not alter tumorigenesis, PTHrP lacking the NLS alone accelerated primary tumor growth by downregulating p27, while PTHrP lacking the NLS and C-terminus repressed tumor growth through p27 induction driven by the tumor suppressor leukemia inhibitory factor receptor (LIFR). Induction of p27 by PTHrP lacking the NLS and C-terminus persisted in bone disseminated cells, but did not prevent metastatic outgrowth, in contrast to the primary tumor site. These data suggest that the PTHrP NLS functions as a tumor suppressor, while the PTHrP C-terminus may act as an oncogenic switch to promote tumor progression through differential regulation of p27 signaling. Introduction Parathyroid hormone-related protein (PTHrP) is a pleiotropic hormone encoded by the PTHLH gene located on chromosome 12, with nine exons and at least three identified promoters [1]. In humans, alternative splicing gives rise to three mature isoforms containing 139, 141, or 173 amino acids, and the first 111 amino acids of the PTHrP sequence are highly conserved among different mammalian species [2]. Regulation of PTHrP is complex and tissue-specific, with the molecule containing numerous cleavage sites and post-translational modifications [1]. The PTHrP polypeptide contains an intracellular *Correspondence: Rachelle W. Johnson rachelle.johnson@vumc.org 1 Graduate Program in Cancer Biology, Vanderbilt University, Nashville, TN, USA 2 Vanderbilt Center for Bone Biology, Vanderbilt University Medical Center, Nashville, TN, USA 3 Department of Medicine, Vanderbilt University Medical Center, Nashville, TN, USA 4 Meharry Medical College, Nashville, TN, USA 5 Department of Cell and Developmental Biology and Program in Developmental Biology, Vanderbilt University, Nashville, TN, USA 6 Vanderbilt University, Nashville, TN 37232, USA Marian University College of Osteopathic Medicine, Indianapolis, IN, USA Bone and Muscle Research Group, Marian University, Indianapolis, IN, USA 9 Academic Affairs, Marian University, Indianapolis, IN, USA 10 Indiana Center for Musculoskeletal Health, Indiana University School of Medicine, Indianapolis, IN, USA 11 Bone Cell Biology and Disease Unit, St. Vincents Institute of Medical Research, Fitzroy, VIC, Australia 12 Department of Medicine, The University of Melbourne, St. Vincents Hospital, Fitzroy, VIC, Australia 7 8 The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the articles Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the articles Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Edwards et al. Breast Cancer Research (2024) 26:34 trafficking and secretion signal, a domain that controls binding to and activation of the classical parathyroid hormone type 1 receptor (PTH1R), and a mid-molecule domain that regulates placental calcium transport. Additionally, the molecule possesses a domain historically termed the nuclear localization sequence (NLS) from amino acids 6794 which regulates nuclear import based on studies carried out in chondrocytes [3], and a carboxy-terminal (C-terminal) domain (beginning at residue 107), to which a number of biological activities have been ascribed [4, 5]. Beyond its well-characterized endocrine and paracrine roles in inducing hypercalcemia of malignancy [6, 7] and tumor-induced bone disease [811], PTHrP regulates the growth of numerous tissues through its intracrine (intracellular) effects on cell survival, proliferation, apoptosis, invasion, and migration, which can occur independent of PTHrP:PTH1R binding on the cell surface [1215]. PTHrP acting through its classical NLS (67-94aa) alters proliferation in peripheral tissues including vascular smooth muscle [1618], where PTHrP also has a smooth muscle relaxing effect [19, 20]. Though less well studied, PTHrP also plays an important role in tumor development. In patients, PTHrP is detectable in most primary breast tumors [11] and serum PTHrP levels are elevated in the majority of patients with hypercalcemia due to breast cancer bone metastases [21, 22]. However, studies have not identified a direct association between elevated serum PTHrP levels in patients and enhanced primary breast tumor growth. The role of PTHrP in primary breast cancer progression remains highly controversial. Some clinical studies demonstrate that PTHrP expression in the primary tumor correlates with improved patient survival and formation of fewer bone metastases [23, 24], while others report that PTHrP is associated with worse patient outcomes [11, 25, 26]. Conflicting data from pre-clinical studies have further confounded the field; genetically similar mouse models that spontaneously form mammary carcinomas have produced directly conflicting results suggesting that PTHrP can inhibit [27] or promote breast tumorigenesis [28]. Thus, the prognostic role for PTHrP in primary breast tumor progression remains largely unclear. In contrast to its uncertain role in the primary tumor, PTHrP has a well-defined deleterious effect on patient outcomes in later stages of disease progression, where its expression drives bone colonization and metastatic tumor growth [11, 26, 29]. Bone disseminated breast cancer cells secrete osteolytic factors like PTHrP, which induces receptor activator of nuclear factor-B ligand (RANKL)-dependent osteoclastogenesis through PTH1R activation on osteoblasts [30]. In human MCF7 breast cancer cells, which normally lie dormant in bone [9, 3133], overexpression of PTHrP (1-139aa) reprograms Page 2 of 17 the cells to become highly osteolytic and dramatically increases bone tumor burden in vivo [9]. Our studies suggest that this potentially occurs through PTHrPmediated suppression of the breast tumor suppressor leukemia inhibitory factor receptor (LIFR) [32, 34, 35] and other pro-dormancy factors [30, 3336]. Our group, and others, have reported evidence that PTHrP can regulate breast tumor progression independent of paracrine or autocrine activation of PTH1R or downstream canonical cAMP signaling [37, 38]. This suggests that PTHrP acts in an intracrine manner to influence breast tumor cell behavior. In support of this, PTHrP (38-94aa) containing the calcium transport region and NLS has been shown to bind to chromatin [39], and full-length secreted PTHrP (-36-139aa) has been shown to localize to the LIFR proximal promoter [40]. In this study, we sought to determine how the intracrine activity of the PTHrP NLS (67-94aa) regulates breast tumor growth and how this effect may be co-regulated by the C-terminal region, since a role for these domains had not been examined in breast cancer cells. In vitro expression of endogenous PTHrP is quite low [32] and there are no reliable antibodies to detect its endogenous isoforms or biological domains. Thus, we rely on an engineered system of expressing truncated mutant proteins with deletion of the PTHrP NLS and C-terminal domains. Our findings begin to provide insight into some of the conflicting preclinical data in the literature, which may provide a framework for targeting PTHrP and its downstream signaling mediators in breast cancer. Results Human breast cancer cells generated to express full-length PTHrP or truncated peptides To determine how PTHrP and its biological domains regulate breast tumor progression, we generated MCF7 human breast cancer cell lines that stably express different domains of the PTHrP molecule (collectively referred to herein as PTHrP mutant cell lines). The plasmids express full-length secreted PTHrP (termed FLSEC, -36-139aa), or truncated forms lacking the classical NLS alone (termed DNLS, -36-6795-139aa) or NLS and C-terminal domain (termed DNLS + CTERM, -36-67aa) with a C-terminal HA tag that is absent in the MSCV control (Fig. 1A). We were unable to generate a mutant with deletion of the secretion signal since these cells do not survive in vitro. We validated plasmid expression at the protein level using an anti-HA antibody and at the mRNA level with qPCR primers targeted to amplify different regions of the Pthlh gene (Fig. 1B-E). To verify expression of the plasmids and characterize the intracellular localization of the PTHrP peptides, we performed immunocytochemical staining for the C-terminal HA tag (Fig. 1F). We confirmed an absence of HA Edwards et al. Breast Cancer Research (2024) 26:34 Page 3 of 17 Fig. 1 Validation of plasmids expressing specific PTHrP domains. (A) Pthlh overexpression construct design and validation in MCF7 cells by (B) western blot for the C-terminal HA-Tag and qPCR for the (C) mid-region, (D) nuclear localization sequence (NLS), and (E) C-terminal domain. MSCV = control, FLSEC = full-length secreted, DNLS = NLS deleted, DNLS + CTERM = NLS and C-terminal domain deleted. Predicted molecular weights: FLSEC PTHrP (-36139aa) = 21.2kD, DNLS PTHrP (-36-67aa)(95-139aa) = 18kD, DNLS + CTERM PTHrP (-36-67aa) = 12.8kD. GAPDH = loading control. (F) Immunocytochemical staining for HA-Tag (green) and DAPI (blue). All panels = 100X and scale bars = 25 m. (G) Secreted PTHrP (1-34aa) levels measured by ELISA from conditioned media of cells described in (A). (B-E & G) n = 3 independent biological replicates. Graphs represent mean SEM. (C) **p < 0.001 vs. MSCV or *p < 0.05 vs. FLSEC by one-way ANOVA with multiple comparisons. (D) **p < 0.001 vs. FLSEC by one-way ANOVA with multiple comparisons. (E) **p < 0.001 vs. DNLS by one-way ANOVA with multiple comparisons expression and fluorescence staining in the MSCV control cells as these plasmids do not contain a C-terminal HA tag. Full-length secreted PTHrP localized to both the nucleus and cytoplasm. Deletion of the NLS alone or NLS and C-terminal domain did not preclude nuclear entry as each PTHrP mutant protein was present in the nucleus as well as cytoplasm (Fig. 1F & Supplementary Fig. 1), regardless of whether they expressed the classical NLS. Therefore, these truncated PTHrP peptides likely gained entry into the nucleus independent of this recognized NLS. While we cannot modulate relative amounts of the PTHrP peptides as it is not possible to accurately Edwards et al. Breast Cancer Research (2024) 26:34 engineer our model system in this manner, we observed no statistically significant difference in PTHrP levels secreted by the PTHrP mutant cell lines compared to controls as measured by an enzyme-linked immunosorbent assay (ELISA) for PTHrP (1-34aa) (Fig. 1G). Thus, altering expression of the NLS or the C-terminal domain does not affect PTHrP secretion by MCF7 cells. Additionally, differences in phenotypes between the PTHrP mutant cells are likely not due to paracrine effects of secreted PTHrP since we and others have previously shown that PTHrP does not activate PTH1R or downstream cAMP signaling in breast cancer cells [37, 38]. The PTHrP NLS and C-terminal domain oppositely regulate breast tumor progression Next, we sought to determine how PTHrP and its biological domains regulate primary breast tumor growth in vivo. Overexpression of full-length PTHrP (-36-139aa) did not significantly alter time to tumor palpation or tumor size compared with controls (Fig. 2A-C and A: p = 0.0497 Log-rank, p = 0.0012 Gehan-Breslow-Wilcoxen; 2B: ANOVA p < 0.0001). Strikingly, deletion of the PTHrP NLS alone resulted in tumors that formed significantly earlier and grew larger than controls, while deletion of both the NLS and C-terminal domains completely reversed this phenotype such that the tumors grew significantly slower and smaller (Fig. 2A-C). To confirm that the PTHrP mutant plasmids were still expressed in vivo, we performed immunofluorescence staining of the primary tumors for the C-terminal HA tag, which was appropriately present in all tumors except the MSCV group, since the MSCV control plasmid does not contain an HA tag (Supplementary Fig. 2). We next assessed whether the changes in tumor size were due to increased proliferation, or reduced cell death. Deletion of the PTHrP NLS alone significantly increased the percentage of Ki67 + positive cells (Fig. 2D) and mitoses (Fig. 2E) in the primary tumors while deletion of both the NLS and C-terminal domain resulted in significantly decreased mitoses (Fig. 2E). There was no difference in cleaved PARP staining in any of the PTHrP mutant cell lines compared to MSCV controls (Fig. 2F). Collectively, these data suggest that the PTHrP NLS regulates breast tumor growth by increasing tumor cell proliferation without impacting apoptosis, but this function is abolished when the PTHrP C-terminus is deleted. p27 is differentially regulated by the PTHrP NLS and C-terminal domains in breast cancer To better understand the in vivo phenotype and mechanism by which the PTHrP NLS and C-terminal domains differentially regulate breast cancer cell proliferation, we performed RNA sequencing on the PTHrP mutant cell lines. We identified several hundred significantly Page 4 of 17 altered genes ( log2 fold change 1 or log2 fold change 1, p < 0.05) that were differentially expressed across the PTHrP mutants (Fig. 3A, Supplementary Data 1). Gene Set Enrichment Analysis (GSEA) of these data revealed that in cells lacking the PTHrP NLS, there was a significant enrichment for genes that are upregulated in MCF7 cells overexpressing the oncoprotein and cell cycle promoter, cyclin D1 (Fig. 3B), indicating that the PTHrP NLS modulates the expression of cell cycle regulators to alter proliferation in MCF7 breast cancer cells. Furthermore, cells expressing PTHrP lacking both the NLS and C-terminal domain were positively enriched for genes involved in the p53 pathway (NES = 1.53, FDR = 0.053). We also examined enriched cancer Hallmark pathways, which revealed an increase in additional cell cycle-related pathways, including G2M Checkpoint and Mitotic Spindle genes (Fig. 3C&D). Based on these RNA sequencing data which pointed to differences in genes encoding cell cycle regulatory proteins, and since p21 and p27 are known to be regulated downstream of PTHrP in other cell types [1618], we investigated these cell cycle factors as a mechanism by which the PTHrP NLS and C-terminal domain oppositely influence breast tumor growth. Immunocytochemical staining revealed that while overexpression of full-length PTHrP (-36-139aa) did not alter p27 levels (Fig. 3E), p27 expression was significantly lower with deletion of the NLS alone compared to control cells. Furthermore, expression of p27 was significantly increased with deletion of both the NLS and C-terminal domain, exceeding levels in both MSCV controls and NLS-alone deleted cells (Fig. 3E). Immunofluorescent staining of the primary breast tumors similarly revealed no change in p27 with overexpression of the full-length PTHrP molecule, but p27 protein levels were significantly decreased with deletion of the NLS alone compared to controls, and oppositely increased with deletion of both the NLS and C-terminal domain (Fig. 3F). Interestingly, in vivo p27 protein levels still remained lower than controls with deletion of both domains (Fig. 3F). When we assessed p21 protein expression, we found inconsistent staining patterns between in vitro cultured cells and in vivo tumor sections; however, we did see a modest increase in p21 staining in tumors expressing full-length secreted PTHrP, suggesting p21 may be regulated downsteam of the intact PTHrP molecule in the context of the tumor microenvironment (Supplementary Fig. 3A&B). Together, these in vitro and in vivo findings suggest that p27 is oppositely regulated by the PTHrP NLS and C-terminal domain in breast cancer, with much lower levels in fast-growing tumors. The difference in p27 expression may therefore contribute to the differential proliferation and breast tumor growth effects observed in vivo. Edwards et al. Breast Cancer Research (2024) 26:34 Page 5 of 17 Fig. 2 Deletion of the PTHrP NLS alters breast cancer cell proliferation and primary tumor growth. (A) Time to tumor palpation, (B) tumor volume over time by digital caliper measurement and (C) final tumor weight in mice inoculated with MSCV, FLSEC, DNLS, or DNLS + CTERM cells into the mammary fat pad. n = 710 mice/group. (D) Ki67 staining and quantification from tumors in (A-C). (E) Quantification of mitoses (# mitotic figures/total cells in 40X field) by DAPI staining from tumors in (A-C). (F) Cleaved PARP staining and quantification from tumors in (A-C). All panels = 40X and scale bar = 50 m. (A) *p < 0.05 vs. MSCV by one-way ANOVA with multiple comparisons or #p < 0.05 DNLS vs. DNLS + CTERM by unpaired t-test. (B) ****p < 0.0001 vs. MSCV by one-way ANOVA with multiple comparisons or **p < 0.01 vs. DNLS by unpaired t-test. (C) **p < 0.01 vs. MSCV by one-way ANOVA with multiple comparisons or ***p < 0.001 vs. DNLS by unpaired t-test. (D) **p < 0.01 vs. MSCV by one-way ANOVA with multiple comparisons. (E) *p < 0.05 vs. MSCV by one-way ANOVA with multiple comparisons or *p < 0.05 vs. DNLS by unpaired t-test. Graphs represent mean SEM PTHrP regulates downstream LIFR signaling to alter p27 expression in vitro We previously demonstrated that PTHrP localizes to the proximal promoter region [40] and downregulates breast cancer cell expression of leukemia inhibitory factor receptor (LIFR) [32], which is a known breast tumor dormancy regulator in bone [32, 36], breast tumor suppressor, and lung metastasis suppressor [34, 35]. The downstream signaling mechanisms by which LIFR regulates breast tumor growth remain incompletely understood. While LIFR is a cell surface receptor, it can also be internalized to the cytoplasm once bound by the LIF ligand [41]. Although overexpression of full-length PTHrP (-36-139aa) has been shown to downregulate Edwards et al. Breast Cancer Research Fig. 3 (See legend on next page.) (2024) 26:34 Page 6 of 17 Edwards et al. Breast Cancer Research (2024) 26:34 Page 7 of 17 (See figure on previous page.) Fig. 3 PTHrP lacking the NLS and C-terminal domain regulates proliferation by altering expression of p27. (A) Number of genes identified by RNAseq with log2fold change > 1 and p < 0.05. (B) GSEA plot from DNLS cells showing enrichment of Cyclin D1 gene signature in MCF7 cells. (C) GSEA plot from DNLS + CTERM cells showing enrichment of genes that regulate the G2M checkpoint. (D) Top twenty enriched Hallmark pathways from FLSEC, DNLS, and DNLS + CTERM cells. (E) Immunocytochemical staining and quantification of p27 in MSCV, FLSEC, DNLS, or DNLS + CTERM cells. n = 3 independent biological replicates. All panels = 40X, scale bar = 25 m. (F) Immunofluorescence staining and quantification for p27 in primary tumors from mice inoculated with MSCV, FLSEC, DNLS, or DNLS + CTERM cells. All panels = 40X, scale bar = 50 m. (E) **p < 0.01 or ****p < 0.0001 vs. MSCV by one-way ANOVA with multiple comparisons or ****p < 0.0001 vs. DNLS by unpaired t-test. (F) *p < 0.05 or ***p < 0.001 vs. MSCV by one-way ANOVA with multiple comparisons or **p < 0.01 vs. DNLS by unpaired t-test. Graphs represent mean SEM LIFR in vitro [32, 36, 38], we observed no difference in LIFR protein expression in vivo with overexpression of the full-length PTHrP molecule (Fig. 4A). Deletion of the PTHrP NLS alone modestly suppressed LIFR levels compared to MSCV controls while deletion of both the NLS and C-terminal domain significantly increased expression of LIFR compared to tumors lacking the NLS alone, which restored levels close to that of the control tumors (Fig. 4A). This pattern of increased LIFR expression with deletion of the PTHrP NLS and C-terminal domain (compared to NLS alone deletion) mirrored the previously observed trend in tumor p27 expression. Thus, we hypothesized that PTHrP may regulate tumor cell proliferation through p27 signaling downstream of LIFR, resulting in altered breast tumor cell proliferation. To investigate this further, we treated the PTHrP mutant cells with a commercially available LIFR inhibitor (EC359) that blocks receptor/ligand interactions. Effective LIFR inhibition was confirmed by decreased phosphorylation of the downstream LIFR signaling factor, pERK (Fig. 4B & D). We did not observe changes in cell cycle phases with LIFR inhibitor treatment of the PTHrP mutant cells (Supplementary Fig. 4A). In the vehicle treated group, p27 remained significantly higher in cells expressing PTHrP lacking the NLS and C-terminal domain compared to those lacking the NLS alone (Fig. 4C). After 24 h of low dose LIFR inhibitor treatment (50nM), this difference was no longer significant (Fig. 4B & C). High dose treatment of LIFR inhibitor (100nM) for 24 h completely reversed the induction of p27 in cells lacking the PTHrP NLS and C-terminal domain such that p27 expression was significantly lower than even control MSCV cells (Fig. 4B & C). Together, these data suggest that PTHrP may induce p27 through a LIFR-dependent mechanism. Treatment of the PTHrP mutant cell lines with the LIFR inhibitor for 1 or 6 h did not elicit the same effect on p27 as the 24-h treatments, such that there was no change in the pattern of p27 protein levels compared with vehicle treated cells (Supplementary Fig. 4B-E). This lack of effect with shorter treatments suggests that p27 is likely an indirect downstream target of LIFR. LIFR is a known dormancy regulator in breast tumor cells in the primary [32, 36, 38] and bone metastatic sites [32]. LIFR signaling activates multiple downstream signaling pathways in breast cancer, including ERK [42]. Since a high p38/ERK signaling ratio promotes tumor dormancy [43, 44], we also analyzed phosphorylated p38 levels in the PTHrP mutant cells, with and without LIFR inhibition. While phosphorylated p38 and the p38/ERK ratio were unchanged in the untreated cells expressing full-length or NLS alone-deleted PTHrP, both p38 and the p38/ERK ratio increased in cells expressing PTHrP lacking the NLS and C-terminal domain, compared to controls (Fig. 4D-F). This suggests that PTHrP lacking the NLS and C-terminal domain preferentially activates p38 signaling, which may induce a more quiescent phenotype. This is consistent with the significantly reduced primary tumor growth (Fig. 2A-C, DNLS + CTERM group). Interestingly, there was a significant increase in phosphorylated p38 and the p38/ERK ratio in the LIFR inhibitor treated cells compared to vehicle treated cells (Fig. 4F). This suggests that the LIFR inibitor may preferentially decrease ERK signaling, which in turn increases p38 activity. Loss of the PTHrP NLS enhances bone metastasis formation despite persistently elevated p27 expression Given the well-established role of PTHrP in promoting metastasis formation [811], we investigated how the NLS and C-terminal domain alter signaling and behavior of bone-disseminated tumor cells using a mouse model of bone colonization in which the PTHrP mutant tumor cells were inoculated through the left cardiac ventricle. We specifically examined whether elevated p27 expression is sustained in bone-disseminated breast tumor cells that express PTHrP lacking the NLS and C-terminal domain and if this alters proliferation, as in the primary tumor. Four weeks post-intracardiac inoculation, qPCR was performed on homogenized femora for human CDKN1B (gene name for p27), and normalized to ACTB (human tumor housekeeping gene) and Hmbs (mouse housekeeping gene) to quantify p27 specifically in bonedisseminated human tumor cells. CDKN1B was significantly higher in the homogenized femora from mice with bone-disseminated tumor cells that expressed PTHrP lacking the NLS and C-terminal domain only (Fig. 5A), confirming that even in the distant metastatic site, the truncated form of PTHrP induces more p27 in tumor cells than other PTHrP peptides. We observed the same trend in p27 expression in the primary tumor. Surprisingly, although p27 levels were higher in the homogenized femora of mice inoculated with tumor Edwards et al. Breast Cancer Research (2024) 26:34 Page 8 of 17 Fig. 4 PTHrP differentially regulates p27 through LIFR in breast cancer cells. (A) Immunofluorescence staining and quantification of LIFR in primary tumors from mice inoculated with MSCV, FLSEC, DNLS, or DNLS + CTERM cells. All panels = 40X and scale bars = 50 m. (B) Western blot analysis of p27, pERK, ERK, p-p38, p38 and tubulin (loading control) protein levels in MSCV, FLSEC, DNLS, or DNLS + CTERM cells treated with vehicle (DMSO) or LIFR inhibitor (EC359, 50nM or 100nM) for 24 h. Densitometry for western blot analysis of (C) p27, (D) pERK/ERK and (E) p-p38/p38 described in (B). (A) **p < 0.01 vs. DNLS by unpaired t-test. (C) *p < 0.05 vs. DNLS by unpaired t-test or *p < 0.05 vs. MSCV by one-way ANOVA with multiple comparisons. (D & E) *p < 0.05 vs. MSCV by one-way ANOVA with multiple comparisons or *p < 0.05, **p < 0.01, ***p < 0.001 versus vehicle by two-way ANOVA. Graphs represent mean SEM cells that express PTHrP lacking the NLS and C-terminal domain, there was significantly elevated osteolytic bone destruction (Fig. 5B-D) and tumor burden (Fig. 5E) in the contralateral limb, as measured by flow cytometric analysis of CD298 + tumor cells, a validated marker for human tumor cells in the bone marrow [45]. The level of metastatic tumor growth and bone destruction was similar in mice inoculated with tumor cells expressing PTHrP either lacking the NLS alone or the NLS and C-terminal domain. This was in striking contrast to the primary tumor site where these cell lines expressing truncated forms of PTHrP elicited opposite effects on breast tumor growth (Fig. 2A-C). Thus, when the NLS and C-terminal domains are deleted, PTHrP induction of p27 persists Edwards et al. Breast Cancer Research (2024) 26:34 Page 9 of 17 Fig. 5 Truncated. PTHrP induces CDKN1B in the bone metastatic site, but enhances osteolysis and tumor burden. (A) qPCR analysis for CDKN1B (p27) normalized to ACTB as a marker of total tumor burden in the bone marrow of mice inoculated with MSCV, FLSEC, DNLS, or DNLS + CTERM cells via intracardiac injection. n = 810 mice/group. (B-D) Total osteolytic lesion area and lesion number (per mouse) based on radiographic analyses for mice described in A. White arrows indicate osteolytic lesions. (E) Flow cytometric quantitation of percent CD298 + tumor cells in the bone marrow of mice described in A. (F) qPCR analysis for RANKL/OPG (Tnfsf11 / Tnfrsf11b) in whole homogenized femurs from mice described in (A). n = 810 mice/group. *p < 0.05, **p < 0.01, or ****p < 0.0001 vs. MSCV by one-way with multiple comparisons. Graphs represent mean SEM in the bone metastatic site. However, in contrast to the primary tumor, induction of p27 downstream of PTHrP in disseminated tumor cells is not sufficient to prevent colonization of the bone and metastatic outgrowth, both of which are elevated by truncated PTHrP peptides lacking the NLS. To determine whether the increase in tumor burden was due to increased osteoclast-mediated bone resorption, we assessed the RANKL/OPG ratio in whole, homogenized femurs across all groups as a marker of osteoclasts. Surprisingly, we only observed a significant increase in RANKL/OPG when the PTHrP NLS domain was deleted, and not in the NLS + C-terminal deleted group. These data suggest that loss of the PTHrP NLS stimulates osteoclast-mediated bone resorption, but loss of the PTHrP NLS and C-terminus does not. We also examined liver histological sections for metastatic tumor burden, but there were no lesions observed in any of the groups. Furthemore, in vitro we observed no difference in migratory potential of cells expressing full-length PTHrP or its truncated forms versus control cells (Supplementary Fig. 5). Together, these data suggest the PTHrP NLS and C-terminal domains may selectively Edwards et al. Breast Cancer Research (2024) 26:34 enhance the ability of breast cancer cells to colonize, survive and proliferate specifically in the bone rather than broadly affecting their ability to migrate from the primary tumor and disseminate to other organs. Discussion PTHrP is a critical driver of tumor-induced bone disease and an important regulator of breast tumorigenesis, cancer progression, and tumor dormancy [28, 32, 46, 47]. Here we investigated the intracellular actions of PTHrP through its NLS and C-terminal domain in breast cancer progression. An important finding is that deletion of the classical PTHrP NLS (67-94aa) does not preclude entry of PTHrP into the nucleus. This indicates that the truncated PTHrP peptides can translocate into the nucleus independent of this recognized NLS. Indeed, one study has reported that PTHrP (1-141) can be endocytosed and translocated into the nucleus via a non-PTH1R cell surface receptor [48], though the mechanism has not been fully elucidated. We are actively investigating alternative mechanisms by which PTHrP enters the nucleus when the classical NLS is deleted. These findings indicate that our study outcomes are likely due to differences in the binding partners or direct interactions of truncated PTHrP with other molecules, rather than the subcellular localization of the truncated peptides. We are also further investigating how the intracellular location alters binding partners of truncated PTHrP peptides to regulate downstream breast cancer cell signaling. Our data demonstrate that the biological domains of PTHrP have distinct functions in breast cancer. These findings are consistent with studies from the skeletal field, which ascribe multiple biological functions to PTHrP domains, particularly through regions outside of the PTH1R-binding domain. Indeed, a knock-in mouse model (PthrpD/D) lacking the midregion, NLS, and C terminal domain (67-137aa) revealed that the intracrine actions of PTHrP are crucial for normal skeletal development and the differentiation of osteogenic and hematopoietic precursors [49]. Most PthrpD/D mice exhibit severe skeletal abnormalities, growth retardation, and die within 5 days. Injection with exogenous PTHrP fails to rescue the lethal phenotype providing further evidence that the effects of PTHrP on these physiological processes are primarily mediated by intracrine signaling. Another in vivo study demonstrated that knock-in mice expressing truncated PTHP (1-84aa) display abnormal skeletal growth and early lethality due to decreased cell proliferation, early senescence, and increased apoptosis in multiple tissues [1618]. Together, these studies demonstrate the importance of the PTHrP NLS and C-terminal domain in regulating tissue development via intracrine signaling, and our data now identify distinct Page 10 of 17 functions of these domains in the pathologic setting of breast cancer. While a large body of evidence indicates that PTHrP has deleterious effects during late stages of breast cancer by promoting bone metastasis, tumor-induced osteolysis, and exit from dormancy, PTHrPs role early in disease progression is highly controversial [27, 28, 32, 46, 47]. Prior preclinical studies reported directly conflicting evidence suggesting that PTHrP inhibits primary breast tumorigenesis in some models [27], while promoting tumor growth in others [28]. Our in vivo findings offer interesting insight into the complex role that PTHrP plays in breast tumor progression. Our data indicate that PTHrP lacking its classical NLS sequence dramatically accelerates breast tumor growth and proliferation in the primary tumor site, suggesting that this domain actually functions to suppress breast tumor growth. Surprisingly, this phenotype is completely reversed if breast cancer cells express PTHrP lacking both the NLS and C-terminal domain, suggesting that the C-terminal domain may possess oncogenic activity that opposes the influence of the NLS. Thus, we are actively pursuing studies to determine how expression or deletion of the C-terminus alone impacts breast cancer growth and bone colonization. Importantly, our data shed light on the conflicting preclinical studies suggesting that PTHrP can promote or inhibit breast tumorigenesis. These controversies may be in part due to the presence of different predominant truncated peptides of PTHrP containing the NLS or C-terminal domain. Unfortunately, these forms are not discernible by commercially available amino-terminal antibodies. While studies have not identified the same engineered fragments as in our model presented here, it is feasible that fragments lacking the classical NLS (67-94aa) or the NLS and C-terminal domain (107-139aa) may naturally circulate in pre-clinical mouse tumor models and patients. In fact, the PTHrP sequence has numerous known and putative mono- and multi-basic cleavage sites [4, 50]. Importantly, PTHrP peptides containing the N-terminal domain (1-36aa), mid-regions (38-94aa), (3895aa) and (38-101aa), as well as the C-terminal domain (107-139aa) have been detected in preclinical mouse models [21, 51] and from the plasma and urine of human patients with solid tumors [21, 51]. While very few studies have investigated a role for these and other PTHrP fragments in breast cancer, some limited studies have identified how their expression alters breast tumor cell behavior, breast tumor growth, and patient outcomes. The PTHrP mid-region fragment (38-94aa) containing a portion of the classical NLS is reported to inhibit in vitro proliferation of MDA-MB-231 human breast cancer cells [52] while another fragment from amino acids 87106 reportedly stimulates proliferation in vitro [53]. Edwards et al. Breast Cancer Research (2024) 26:34 In patients with breast cancer, loss of nuclear localized but not cytoplasmic PTHrP in the primary site has been associated with poor clinical outcomes [54]. Another study identified PTHrP (1248) as a predictive biomarker of breast cancer bone metastasis such that levels of the peptide were significantly increased in the plasma of patients with clinical evidence of bone metastases versus patients without [55]. Together, these studies provide further evidence of domain-specific selectivity for how PTHrP and its truncated isoforms function in vitro versus in vivo. While there were no changes in cell cycling observed in vitro, our in vivo studies demonstrate a modest increase in proliferation with deletion of the NLS alone, which persisted in the primary tumor but not bone. These differences in proliferation in vitro versus in vivo may also be attributed to PTHrP-induced signaling changes in the breast cancer cells that alter their interaction with surrounding stromal cells, including recruitment of immune cells into the tumor microenvironment, which vary substantially by tumor site. The present study sheds important light on the biological role for the classical NLS and C-terminal domain in regulating breast tumor growth in vivo. Examination of cleaved PARP in the primary tumor demonstrated no alterations in apoptosis underlying the differences in tumor burden with expression of PTHrP lacking the NLS alone or both the NLS and C-terminal domain. We also examined levels of cleaved caspase-3 to more broadly assess apoptosis. One limitation in our model is that expression of caspase-3 is low at baseline in MCF7 cells, making it difficult to detect further reductions, particularly in cells expressing PTHrP lacking the DNLS. Importantly, the tumors assessed in our study were analyzed at endpoint, but it is possible that more dramatic changes in apoptosis occurred early in tumor progression. Indeed, the majority of tumors expressing PTHrP lacking the NLS and C-terminal domain were small in size and nearly undetectable at endpoint. The ability to measure apoptotic or proliferative markers from all tumors may have demonstrated a greater difference to further explain the alterations in tumor burden. Cyclin dependent kinase inhibitor proteins are regulated downstream of the PTHrP NLS and C-terminal domain in non-breast cancer cell lineages [1618]. Our studies demonstrate that p27 is oppositely regulated by the PTHrP NLS and C-terminal domain in breast cancer and may be an important downstream signaling factor mediating how these domains differentially alter breast tumor growth (Fig. 6). Specifically, the PTHrP C-terminal domain appears to function as an oncogenic molecular switch able to induce proliferation and promote primary breast tumor formation through a partially LIFRdependent mechanism that suppresses p27 expression. Page 11 of 17 It should be noted that there are significant differences in tumor burden and p27 between control tumors and tumors expressing PTHrP that lack the NLS, but a nonsignificant decrease in LIFR (~ 50% reduction). Thus, the data are consistent across our in vivo study, but do not always result in statistically significant changes. This suggests that LIFR is not the only driver of p27 in our model. Future studies utilizing breast cancer cells expressing PTHrP with deletion of the C-terminal domain only will be needed to confirm this. Interestingly, although CDKN1B (gene name for p27) remained elevated by the bone-disseminated tumor cells expressing PTHrP lacking the NLS and C-terminal domain, the cells readily colonized the bone marrow. We thought this may be due to an increase in osteoclast-mediated bone resorption, which we assessed by measuring RANKL/OPG levels in whole, homogenized femora. We were surprised that RANKL/OPG was only elevated when the PTHrP NLS was deleted, and not when the NLS and C-termainal domain were deleted, since both groups had similar levels of bone destruction and bone metastatic tumor burden. This finding suggests that the mechanism of tumor outgrowth caused by the PTHrP fragments is likely distinct, and that the osteoclast-mediated osteolysis must have occurred early in disease progression in the tumors lacking the PTHrP NLS and C-terminus, since measurements were assessed at endpoint. Follow-up studies to identify the distinct mechanisms of tumor outgrowth in bone that are caused by each PTHrP fragment are underway. In our studies, pharmacologic LIFR inhibition revealed an unexpected trend whereby breast cancer cells treated with the inhibitor had significantly elevated phosphorylated p38 and a p38/ERK signaling ratio compared to vehicle treated cells, regardless of PTHrP mutant expression. This effect was further elevated when the PTHrP NLS and C-terminal domain were deleted. LIFR is known to activate STAT3, ERK, and AKT signaling, among numerous other signaling pathways in breast cancer [32, 42, 56]. It has been postulated that LIFR signaling promotes tumor dormancy specifically through STAT3 activation [32]; however, the oncogenic ERK and AKT pathways can still be activated by LIFR-binding cytokines [42]. Our data here suggest that the EC359 LIFR inhibitor may preferentially decrease LIFR-mediated ERK signaling, shifting the balance towards p38 activity and suppression of cell proliferation in vitro. Since LIFR activates multiple singaling pathways in breast cancer cells [42], we also sought to analyze alterations in STAT3 and AKT signaling in the presence and absence of LIFR inhibition via western blot analysis; however, activation of these pathways was too low at baseline to quantify discernable changes in pSTAT3 and pAKT. Recently, small molecule inhibitors and neutralizing antibodies targeting LIFR have been investigated as a Edwards et al. Breast Cancer Research (2024) 26:34 Page 12 of 17 Fig. 6 Model of PTHrP domain-specific actions in breast cancer progression and bone colonization. In the primary breast site (top left panel, left of arrows), PTHrP lacking the NLS and C-terminal domain decreases tumor cell proliferation through p27 induction driven by the tumor suppressor leukemia inhibitory factor receptor (LIFR). PTHrP lacking the NLS and C-terminal domain also preferentially induces p38 phosphorylation and signaling to inhibit cell cycling downstream of LIFR activation. In the breast, truncated PTHrP lacking the NLS alone (top left panel, right of arrows) downregulates LIFR expression (denoted by transparent coloring) and prevents induction of p27 expression and activation of p38 signaling (denoted by dashed arrows, dotted outlines and transparent coloring) to drive cell proliferation and tumor growth. In bone disseminated tumor cells (bottom panel), LIFR expression is downregulated and the induction of p27 by PTHrP lacking the NLS and C-terminal domain persists, but is not sufficient to repress metastatic outgrowth (denoted by dashed inhibitor line), in contrast to the primary tumor. In the bone, tumor cells expressing PTHrP peptides lacking the NLS or NLS and Cterminal domain readily proliferate into metastatic tumors. Image created with Biorender.com strategy to inhibit breast tumor growth and metastasis in preclinical studies [57, 58]. Although anti-LIFR agents do show evidence of effectively limiting primary breast tumor growth, caution should still be exercised in their use as a breast cancer therapy since inhibiting LIFR signaling could inadvertently increase metastatic outgrowth in bone where the LIFR:STAT3 pathway suppresses proliferation of disseminated breast tumor cells [32, 5961]. It will therefore be important to define the downstream pathways that are disrupted by individual LIFR antagonists. Furthermore, it is still unclear how the PTHrP NLS and C-terminal domains may differentially regulate other downstream LIFR signaling pathways. Concluding remarks In summary, these data reveal important insights into how the PTHrP NLS and C-terminal domain divergently control breast cancer progression through p27 signaling in the primary tumor and bone metastatic site. As a potent regulator of breast tumor growth and distant metastatic progression, PTHrP has the potential to be leveraged as a therapeutic target for the treatment of breast cancer at multiple stages of disease progression and possibly for the prevention of bone metastasis formation. However, it is critical that this work be approached with attention to the PTHrP peptides present and their ability to differentially activate downstream signaling pathways. Edwards et al. Breast Cancer Research (2024) 26:34 Materials and methods Cell culture and reagents Cells PTHrP mutant cell lines were established in the laboratory of one of us (TJM) at St. Vincents Institute of Medical Research, as previously described [61]. Briefly, the following constructs were synthesized by Integrated DNA Technologies (IDT) (Coralville, IA, USA): Pthlh(36-139), Pthlh (1-139), Pthlh(-36-67), Pthlh(-36-139). Xho1/ EcoR1 enzyme digestion and ligation was performed to clone the constructs into the murine stem cell virus (MSCV)-zeo plasmid. Each plasmid except for the MSCV control was tagged with a human influenza hemagglutinin (HA) epitope at the C-terminal end. DNA sequencing was performed by the Australian Genome Research Facility. Phoenix cells were then transfected with the mutant plasmids and used to infect MCF7 cells which were placed under antibiotic selection with Zeocin to establish stable lines. The resulting PTHrP mutant cells were cultured in DMEM supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin (P/S). All cell lines were regularly tested for mycoplasma contamination. Proliferation assays Cells were plated at 1 106 cells per 10cm2 plate and allowed to adhere for 46 h. Adherent cells were then trypsinized and mixed with 0.4% trypan blue solution. Viable cells were determined based on dye exclusion and counted using a TC20 Automated Cell Counter (BioRad). Proliferation of PTHrP mutant cells was monitored daily for four days by trypsinizing and counting viable cells. LIFR inhibitor treatment Cells were plated at 1 106 cells/ 10cm2 plate and allowed to adhere overnight. The following day, cells were treated with EC359, a leukemia inhibitory factor receptor (LIFR) inhibitor (50nM or 100nM; MedChemExpress; Catalog No. HY-1,201,420) or vehicle (0.1% dimethyl sulfoxide, DMSO) for 1, 6, or 24 h in full-serum media. RNA extraction and real-time qPCR RNA was extracted from cells using TRIzol (ThermoFisher) and prepared for real-time qPCR analysis as previously described [32]. Human primers for b2M [32] and CDKN1B (p27) were previously published. The following primers were designed using PrimerBlast (NCBI) against the human genome and validated by dissociation: ACTB (F- CATGTACGTTGCTATCCAGGC), R- CTCCTTAA TGTCACGCACGAT). Mouse primers for HMBS were previously published [32]. The following primers were designed using PrimerBlast (NCBI) against the mouse genome (Mus musculus) and validated by dissociation: Page 13 of 17 PTHrP mid-region (F- CATCAGCTACTGCATGACA AGG, R- GGTGGTTTTTGGTGTTGGGTG), PTHrP NLS (F- AACAGCCACTCAAGACACCC, R- GACCGA GTCCTTCGCTTCTT), PTHrP C-terminal region (F- A AAAGAAGCGAAGGACTCGG, R- GCGTCCTTAAGC TGGGCT). Western blotting Cultured cells were rinsed twice with cold 1X PBS and harvested in RIPA lysis buffer (Sigma) containing protease and phosphatase inhibitors (Roche). Protein lysate (20g) was loaded onto an SDS-PAGE gel under reducing conditions and transferred to nitrocellulose membranes. Membranes were probed with antibodies against HA-Tag (Cell Signaling, C29F4, Catalog No. 37T4S, 1:1000), LIFR (Santa Cruz, C-19, Catalog No. sc-659, 1:1000), p21Waf1/ Cip1(Cell Signaling, Catalog No. 2947 S, 1:1000), p27 Kip1 (Cell Signaling, Catalog No. 3686 S, 1:1000), phospho-p38 MAPK (Thr180/Tyr182) (Cell Signaling, Catalog No. 4511, 1:1000), p38 MAPK (Cell Signaling, Catalog No. 8690, 1:1000), phospho-ERK1/2 Thr202/Tyr204 (Cell Signaling, Catalog No. 9101, 1:1000), ERK1/2 (Cell Signaling, catalog number 9102, 1:1000), Calnexin (AbCam, Catalog No. ab22595-100UG, 1:900), GAPDH (Cell Signaling 14C10, Catalog No. 2118 S, 1:5000), HDAC2 (Cell Signaling, D6S5P, 1:1000), -tubulin (Antibody & Protein Resource at Vanderbilt University, Catalog No. VAPRTUB, 1:5000), or Vinculin (Millipore, Catalog No. AB6039, 1:1000). Nuclear and cytoplasmic extraction Nuclear and cytoplasmic extracts were obtained from cultured PTHrP mutant cells using the NE-PER Nuclear and Cytoplasmic Extraction Reagents Kit (Thermo Scientific, Catalog No. 78,835) according to the manufacturers instructions. Briefly, 5 106 cells were plated in full serum DMEM and allowed to adhere overnight. The following day, adherent cells were trypsinized and centrifuged at 500 x g for 5 min, and the pellet was suspended in PBS. Cells were then transferred to a new microcentrifuge tube and centrifuged at 500 x g for 3 min. Supernatant was discarded and 500 l of ice-cold CER I with 5 l of protease inhibitor was added to the cell pellet and vortexed. The cell suspension was incubated on ice for 10 min. Ice-cold CER II (27.5 l) was then added to the tube, vortexed, and incubated on ice for 1 min. Next, the sample was vortexed and centrifuged at 16,000 x g for 5 min. The supernatant (cytoplasmic extract) was immediately transferred to a clean pre-chilled tube and stored at -80oC. The cell pellet was suspended in 250 l of icecold NER, vortexed for 15 s, and placed on ice. Vortexing was repeated every 10 min for a total of 40 min. The tube was then centrifuged at 16,000 x g for 10 min. Finally, the Edwards et al. Breast Cancer Research (2024) 26:34 supernatant (nuclear extract) was transferred to a clean pre-chilled tube and stored at -80oC. Immunocytochemistry For analysis of HA-tagged PTHrP peptides, cells were seeded onto a 4-well culture slide at 6 105 cells/ well and allowed to adhere overnight. The following day cells were washed twice with 1x PBS and fixed with 10% formalin for 15 min. Cells were then washed three times with 1X PBS for 5 minutes each, permeabilized in 0.25% TritonX in 1X PBS for 10 min and washed twice with 1X PBS for 5 minutes each. Next cells were blocked in a 3% mix of donkey horse serum (DHS)/ bovine serum albumin (BSA) for 1 h at room temperature, washed twice with 1X PBS for 5 minutes each and finally incubated with HATag antibody (Cell Signaling, C29F4, Catalog No. 37T4S, 1:500) diluted in DHS/ BSA mix for 1 h at room temperature. Afterwards, cells were washed three times with 1X PBS for 5 minutes each and incubated in goat anti-rabbit IgG (H + L) Alexa Fluor 488 secondary antibody (Thermo Fisher, Catalog No A-11,034, 1:1000) diluted in DHS/ BSA mix in the dark for 1 h at room temperature. Cells were then washed three times with 1X PBS for 5 minutes each. Lastly, the chamber was removed from each slide before mounting coverslips with VECTASHIELD HardSet Antifade Mounting Medium with DAPI (Vector Laboratories). Fixed cells were imaged on a laser scanning confocal microscope Nikon A1r based on a TiE motorized Inverted Microscope using a (I) 60X lens, NA 1.4, run by NIS Elements C software with sections imaged in 0.23 m slices or (II) 100X lens, NA 1.49, run by NIS Elements C software with sections imaged in 0.23 m slices. For analysis of p21 and p27, 8 105 cells were seeded onto glass coverslips coated with 5 g/ml human fibronectin (Millipore) 12 h prior. The following day, cells were washed with 1X PBS, fixed with 10% formalin for 15 min, washed three times with 1X PBS for five minutes each and permeabilized with 0.25% Triton-X for 10 min. Afterwards, cells were washed twice with 1X PBS for 5 minutes each and blocked with DHS/ BSA mix for 1 h at room temperature. Cells were then washed twice with 1X PBS for 5 minutes each and incubated in p21Waf1/ Cip1(Cell Signaling, Catalog No. 2947 S, 1:1000) or p27 Kip1 (Cell Signaling, Catalog No. 3686 S, 1:1000) diluted in DHS/BSA mix for 1.5 h at room temperature. Afterwards cells were washed three times with 1X PBS for 5 minutes each and incubated in goat anti-rabbit IgG (H + L) Alexa Fluor 488 secondary antibody (Thermo Fisher, Catalog No A-11,034, 1:1000) diluted in DHS/ BSA mix in the dark at room temperature. Finally, cells were washed three times with 1X PBS for 5 minutes each before mounting on glass slides with VECTASHIELD HardSet Antifade Mounting Medium with DAPI (Vector Laboratories). Images were collected on an Olympus Page 14 of 17 BX41 Microscope equipped with an Olympus DP71 camera using the 40X plain objective. For p21 quantitation in Image J, total nuclei and positive staining cells were counted manually to calculate the percent of positive staining cells. For p27, the fluorescence intensity was quantified using ImageJ with manual cell contouring and measurement of the Raw Integrated Density which was averaged across all cells from 3 separate images. Enzyme-linked immunosorbent assay To prepare conditioned media, PTHrP mutant cells (1 105) were plated in full-serum media in a 24-well plate and allowed to adhere for 24 h. Afterwards, the fullserum media was changed to 600 l of reduced serum media (DMEM + 2% FBS + 1% P/S) and cells were incubated for 24 h. Conditioned cell media was then harvested and centrifuged at 1500 rpm for 10 min at 4 C. The supernatant was treated with protease inhibitor (Sigma, P8340, 1:100) before further analysis. Undiluted conditioned media was added to 96-well ELISA plates to measure secreted PTHrP levels according to the manufacturers protocol (Creative Diagnostics, Catalog No. DEIA2034). For the final analysis, calculated PTHrP concentrations measured by the ELISA were normalized to the total protein concentration (mg/ml) in each sample measured by BCA assay (Thermo Fisher). Cell cycle analysis Cell cycle analysis was performed by seeding 150,000 cells per well into 6-well plates for each cell line. After 24 h, cells were treated with 50nM EC359, 100nM EC359, or DMSO vehicle for 48 h. After 48 h, 150,000 cells were removed from each treatment group and live stained with Hoescht 33342 (AbCam) at a concentration of 10 g/mL for 1 h at 37 C. Stained cells were analyzed on a 4 Laser Fortessa by the Vanderbilt Flow Cytometry Resource Core. Flow cytometer data were analyzed using FlowJo software to gate for G0/1, S, and G2 phases. Each bar represents data from 3 independent experiments. Migration assay Scratch assays were performed by seeding 400,000 cells of each mutant cell line (MSCV, FLSEC, DNLS, and DNLS + CTERM) into one well of a 6-well plate. After 24 h, three scratches were made in each well with a pipette tip. Images were taken at 100x on an inverted microscope at 0 h (immediately after scratch), 24 h, and 48 h. Percent closure was determined via analysis with ImageJ. Each replicate is expressed as an average of three scratches per well. Each data point represents three independent experiments. Edwards et al. Breast Cancer Research (2024) 26:34 Animal studies and imaging Animals Experiments were performed under the regulations of the Animal Welfare Act and the Guide for the Care and Use of Laboratory Animals and approved by the Vanderbilt University Institutional Animal Care and Use Committee (IACUC). For the mammary fat pad study, 17-estradiol pellets (0.36 mg/pellet; Innovative Research of America, Catalog No. SE-121) were subcutaneously implanted into female athymic nude mice 24 h prior to tumor inoculation [61]. The following day, 5 105 tumor cells from each pooled cell line in 20 l PBS + 50% matrigel (Fisher Scientific) were inoculated into the fourth mammary fat pad (n = 10 mice injected per group). Tumor volume was assessed by caliper measurement. Several mice had to be sacrificed early due to estrogen-induced toxicities resulting in MSCV = 8 mice, FLSEC = 7 mice, DNLS = 10 mice, DNLS + CTERM = 9 mice in the final analysis. For the intracardiac inoculation study, 6-week-old female athymic nude mice (Jackson, Catalog No. 7850) were injected with 1 105 tumor cells from each pooled cell line as previously described [63] (n = 810 mice injected per group). The mice were subcutaneously implanted with a slow-release 17-estradiol pellet (0.36 mg/pellet; Innovative Research of America, Catalog No. SE-121) 24 h prior to tumor cell injection [63]. Radiography Radiographic (x-ray) images were obtained as previously described [64]. Briefly, a Faxitron LX-60 (34 kV for 8 s) was used to acquire x-ray images and images were quantified for osteolytic lesion number and area using ImageJ software. Histology Upon sacrifice of the mice, dissected tumors were fixed in 10% formalin for 48 h and stored in 70% ethanol until being paraffin-embedded for further analyses. Tissue sections were deparaffinized by heating the slides to 50 C and placed in xylene for 5 min and then 3 min. Next, slides were soaked in 100%, 95%, and then 75% ethanol for 3 min each. Slides were slowly changed to deionized water and rinsed twice in water. The slides were immersed in 10 mM TRIS (pH 9.0) and 1 mM EDTA heated to 150 C for 20 min. After cooling at room temperature for 20 min, slides were rinsed twice with water and then three times with 1X PBS followed by blocking with 10% BSA in PBS for 2 h. Sections were stained with Ki67 (Thermo Fisher; Catalog No. RM9106S0, 1:500), cleaved PARP (Asp214) (Cell Signaling Technology, Catalog No. 5625T, 1:500), HA-Tag (Cell Signaling, C29F4, Catalog No. 37T4S, 1:1000), p21Waf1/Cip1(Cell Signaling, Catalog No. Page 15 of 17 2947 S, 1:1000), or p27 Kip1 (Cell Signaling, Catalog No. 3686 S, 1:1000) in 3% BSA in PBS overnight at 4 C. The following day, sections were washed three times with 1X PBS and incubated in goat anti-rabbit IgG (H + L) Alexa Fluor 488 secondary antibody (Thermo Fisher, Catalog No A-11,034, 1:1000) in 3% BSA/PBS in the dark at room temperature for 1 h. Finally, sections were washed three times with 1X PBS and coverslips were mounted using VECTASHIELD HardSet Antifade Mounting Medium with DAPI (Vector Laboratories). For LIFR staining, after blocking in 10% BSA for 2 h, slides were incubated in FITC-LIFR (Santa Cruz, Catalog No. sc-515,337, 1:50) in 3% BSA/PBS overnight at 4 C. The following day, sections were washed three times with 1X PBS and coverslips mounted using VECTASHIELD HardSet Antifade Mounting Medium with DAPI (Vector Laboratories). All images except for Ki67 were collected on an Olympus BX41 Microscope equipped with an Olympus DP71 camera using the 40X plain objectives. For LIFR quantitation, 40X images were used and an area measuring 1900 1180 pixels was selected to measure the Raw Integrated Density. The Raw Integrated Density from 3 representative images was averaged for each mouse and these values are reported in the figure. For p21, p27, and cleaved PARP, the quantitation was performed using ImageJ analysis of the 40X images. Positive staining nuclei and total cell counts were determined using color thresholding in ImageJ and the number of positive staining nuclei was divided by the total number of nuclei present to calculate the percent positivity. For Ki67 quantification, fixed samples were imaged on a laser scanning confocal microscope Nikon A1r based on a TiE motorized Inverted Microscope using a 60X lens, NA 1.4, run by NIS Elements C software. Sections were imaged in 0.4 m slices. Positive staining nuclei and cell counts were determined using color thresholding in ImageJ and the number of positive staining nuclei was divided by the total number of nuclei present to calculate percent Ki67 positivity. Flow Cytometry One hindlimb (inclusive of bone marrow and tumor cells) was crushed with a mortar and pestle to obtain the bone marrow. PBS (1mL) was added to the crushed bone marrow and spun down and washed with PBS to remove bone debris. Bone marrow (5 105 cells) was stained in 100L of PBS with LIVE/DEAD Fixable Green Dead Cell Stain Kit @488nm (Thermo Fisher Scientific, Catalog Number L34970, 1:1000) for 15 min on ice at 4 C in the dark. Cells were washed with PBS and resuspended with 100L of 1% BSA in PBS with CD298 antibody (BioLegend, Cat #341,704) for 30 min on ice at 4 C in the dark. Edwards et al. Breast Cancer Research (2024) 26:34 Flow Cytometry Analysis Flow cytometry experiments were performed in the VUMC Flow Cytometry Shared Resource using the 5-laser BD LSRII and 4-laser BD Fortessa LSRII. Data was analyzed using FlowJo software (FlowJo, LLC) where bone marrow samples were gated based on forward scatter and side scatter geometry, and PE-CD298 (+) cells were gated using live cells (LIVE/DEAD-Green negative) as previously validated in tumor-bearing bone marrow samples [45]. MCF7 breast cancer cells were used as a positive control for CD298 stain. Statistics and reproducibility For all experiments, n per group is as indicated by the figure legend and the scatter dot plots indicate the mean of each group and error bars indicate the standard error of the mean. All graphs and statistical analyses were generated using Prism software (Graphpad). Statistical significance for all in vitro and in vivo assays was analyzed using an unpaired t-test, one-way ANOVA with Sidaks multiple comparisons test or two-way ANOVA with multiple comparisons, as indicated in the figure legends. For each analysis p < 0.05 was considered statistically significant, and *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. Supplementary Information The online version contains supplementary material available at https://doi. org/10.1186/s13058-024-01791-z. Supplementary Material 1 Supplementary Material 2 Supplementary Material 3 Author contributions C.E.M. wrote the main manuscript text, performed experiments, analyzed data, and prepared all figures. J.F.K., J.A.S., D.M.G, J.A.J, M.A.H.D., L.A.V.III, K.M.B., T.N.O., J.R.F., B.A.K., H.T.S., and C.J.V. performed experiments and analyzed data included in Figs. 1, 2, 3, 4 and 5. J.W.L. and T.J.M. generated resources and reagents for experiments and provided project input. R.W.J. wrote and edited the manuscript text, prepared figures, analyzed data, and secured funding for the project. All authors reviewed the manuscript. Funding This work was supported by DoD Breakthrough Award W81XWH-22-1-0090 (R.W.J.). This project was also supported by scholarship funds from NIH award P30CA06848 Vanderbilt-Ingram Cancer Center Support Grant and NIGMS T32GM007347. Data availability Data that support the findings of this study are available from the corresponding author upon reasonable request. Declarations Ethical approval Experiments were performed under the regulations of the Animal Welfare Act and the Guide for the Care and Use of Laboratory Animals and approved by the Vanderbilt University Institutional Animal Care and Use Committee (IACUC). Page 16 of 17 Competing interests The authors declare no competing interests. Received: 14 July 2023 / Accepted: 19 February 2024 References 1. Martin TJ, Moseley JM, Gillespie MT. Parathyroid hormone-related protein: biochemistry and molecular biology. Crit Rev Biochem Mol Biol. 1991;26(34):37795. 2. McCauley LK, Martin TJ. Twenty-five years of PTHrP progress: from cancer hormone to multifunctional cytokine. J Bone Min Res. 2012;27(6):12319. 3. Henderson JE, et al. Nucleolar localization of parathyroid hormone-related peptide enhances survival of chondrocytes under conditions that promote apoptotic cell death. Mol Cell Biol. 1995;15(8):406475. 4. Martin TJ. Parathyroid hormone-related protein, its regulation of cartilage and bone development, and role in treating Bone diseases. Physiol Rev. 2016;96(3):83171. 5. Martin TJ, Sims NA et al. Chap. 25 - Paracrine parathyroid hormonerelated protein in bone: physiology and pharmacology, in Principles of Bone Biology (Fourth Edition), J.P. Bilezikian, Editors. 2020, Academic Press. p. 595621. 6. Albright F, Smith PH, Richardson AM. Postmenopausal osteoporosis: its clinical features. JAMA. 1941;116(22):246574. 7. Burtis WJ, et al. Identification of a novel 17,000-dalton parathyroid hormonelike adenylate cyclase-stimulating protein from a tumor associated with humoral hypercalcemia of malignancy. J Biol Chem. 1987;262(15):71516. 8. Powell GJ, et al. Localization of parathyroid hormone-related protein in breast cancer metastases: increased incidence in bone compared with other sites. Cancer Res. 1991;51(11):305961. 9. Thomas RJ, et al. Breast cancer cells interact with osteoblasts to support osteoclast formation. Endocrinology. 1999;140(10):44518. 10. Guise TA, et al. Evidence for a causal role of parathyroid hormone-related protein in the pathogenesis of human breast cancer-mediated osteolysis. J Clin Invest. 1996;98(7):15449. 11. Southby J, et al. Immunohistochemical localization of parathyroid hormonerelated protein in human breast cancer. Cancer Res. 1990;50(23):77106. 12. Dougherty KM, et al. Parathyroid hormone-related protein as a growth regulator of prostate carcinoma. Cancer Res. 1999;59(23):601522. 13. Phan TG, Croucher PI. The dormant cancer cell life cycle. Nat Rev Cancer. 2020;20(7):398411. 14. Shen X, et al. Increased cell survival, migration, invasion, and akt expression in PTHrP-overexpressing LoVo colon cancer cell lines. Regul Pept. 2007;141(13):6172. 15. Tovar Sepulveda VA, Shen X, Falzon M. Intracrine PTHrP protects against serum starvation-induced apoptosis and regulates the cell cycle in MCF-7 breast cancer cells. Endocrinology. 2002;143(2):596606. 16. de Miguel F, et al. The C-Terminal region of PTHrP, in Addition to the Nuclear Localization Signal, is essential for the Intracrine Stimulation of Proliferation in Vascular smooth muscle cells. Endocrinology. 2001;142(9):4096105. 17. Fiaschi-Taesch N, et al. Mutant parathyroid hormone-related protein, devoid of the nuclear localization signal, markedly inhibits arterial smooth muscle cell cycle and neointima formation by coordinate up-regulation of p15Ink4b and p27kip1. Endocrinology. 2009;150(3):142939. 18. Miao D et al. Severe growth retardation and early lethality in mice lacking the nuclear localization sequence and C-terminus of PTH-related protein Proceedings of the National Academy of Sciences, 2008;105(51):20309. 19. Maeda S, et al. Targeted overexpression of parathyroid hormone-related protein (PTHrP) to vascular smooth muscle in transgenic mice lowers blood pressure and alters vascular contractility. Endocrinology. 1999;140(4):181525. 20. Nishikawa N, et al. PTHrP is endogenous relaxant for spontaneous smooth muscle contraction in urinary bladder of female rat. Endocrinology. 2013;154(6):205868. 21. Burtis WJ, et al. Preliminary characterization of circulating amino- and carboxy-terminal fragments of parathyroid hormone-related peptide in humoral hypercalcemia of malignancy. J Clin Endocrinol Metab. 1992;75(4):11104. 22. Grill V, et al. Parathyroid hormone-related protein: elevated levels in both humoral hypercalcemia of malignancy and hypercalcemia complicating metastatic breast cancer. J Clin Endocrinol Metab. 1991;73(6):130915. Edwards et al. Breast Cancer Research (2024) 26:34 23. Henderson MA, et al. Parathyroid hormone-related protein production by breast cancers, Improved Survival, and reduced bone metastases. JNCI: J Natl Cancer Inst. 2001;93(3):2347. 24. Tran TH et al. Loss of Nuclear localized parathyroid hormone-related protein in primary breast Cancer predicts poor clinical outcome and correlates with suppressed Stat5 signaling. 2018;24(24):635566. 25. Linforth R, et al. Coexpression of parathyroid hormone related protein and its receptor in early breast cancer predicts poor patient survival. Clin Cancer Res. 2002;8(10):31727. 26. Yoshida A, et al. Significance of the parathyroid hormone-related protein expression in breast carcinoma. Breast Cancer. 2000;7(3):21520. 27. Fleming, N.I., et al., Parathyroid HormoneRelated Protein Protects against Mammary Tumor Emergence and Is Associated with Monocyte Infiltration in Ductal Carcinoma In situ Cancer Research, 2009;69(18):7473. 28. Li J, et al. PTHrP drives breast tumor initiation, progression, and metastasis in mice and is a potential therapy target. J Clin Investig. 2011;121(12):465569. 29. Bundred NJ, et al. Parathyroid hormone related protein and skeletal morbidity in breast cancer. Eur J Cancer. 1992;28(23):6902. 30. Mundy GR. Mechanisms of bone metastasis. Cancer. 1997;80(8 Suppl):154656. 31. Barkan D, et al. Inhibition of metastatic outgrowth from single dormant tumor cells by targeting the cytoskeleton. Cancer Res. 2008;68(15):624150. 32. Johnson RW, et al. Induction of LIFR confers a dormancy phenotype in breast cancer cells disseminated to the bone marrow. Nat Cell Biol. 2016;18(10):107889. 33. Wang H, et al. The osteogenic niche promotes early-stage bone colonization of disseminated breast cancer cells. Cancer Cell. 2015;27(2):193210. 34. Chen D, et al. LIFR is a breast cancer metastasis suppressor upstream of the Hippo-YAP pathway and a prognostic marker. Nat Med. 2012;18(10):15117. 35. Iorns E, et al. Whole genome in vivo RNAi screening identifies the leukemia inhibitory factor receptor as a novel breast tumor suppressor. Breast Cancer Res Treat. 2012;135(1):7991. 36. Clements ME et al. HDAC inhibitors induce LIFR expression and promote a dormancy phenotype in breast cancer. Oncogene, 2021. 37. Grinman DY, et al. PTHrP induces STAT5 activation, secretory differentiation and accelerates mammary tumor development. Breast Cancer Res. 2022;24(1):30. 38. Johnson RW, et al. Parathyroid hormone-related protein negatively regulates Tumor Cell Dormancy genes in a PTHR1/Cyclic AMP-Independent manner. Front Endocrinol. 2018;9:2411. 39. Sirchia R, et al. Mid-region parathyroid hormone-related protein (PTHrP) binds chromatin of MDA-MB231 breast cancer cells and isolated oligonucleotides in vitro. Breast Cancer Res Treat. 2007;105(1):10516. 40. Edwards CM, et al. HDAC inhibitors stimulate LIFR when it is repressed by hypoxia or PTHrP in breast cancer. J Bone Oncol. 2021;31:100407. 41. Hilton DJ, Nicola NA. Kinetic analyses of the binding of leukemia inhibitory factor to receptor on cells and membranes and in detergent solution. J Biol Chem. 1992;267(15):1023847. 42. Omokehinde T, Jotte A, Johnson RW. gp130 cytokines activate Novel Signaling pathways and alter bone dissemination in ER + breast. Cancer Cells. 2022;37(2):185201. 43. Aguirre-Ghiso JA, et al. ERK(MAPK) activity as a determinant of tumor growth and dormancy; regulation by p38(SAPK). Cancer Res. 2003;63(7):168495. 44. Bragado P, et al. TGF-beta2 dictates disseminated tumour cell fate in target organs through TGF-beta-RIII and p38alpha/beta signalling. Nat Cell Biol. 2013;15(11):135161. 45. Sowder ME, Johnson RW. Enrichment and detection of bone disseminated tumor cells in models of low tumor burden. Sci Rep. 2018;8(1):14299. Page 17 of 17 46. Edwards CM, Johnson RW. From Good to Bad: The Opposing Effects of PTHrP on Tumor Growth, Dormancy, and Metastasis Throughout Cancer Progression. 2021;11(728). 47. Martin TJ, Moseley JM. Mechanisms in the skeletal complications of breast cancer. Endocr Relat Cancer. 2000;7(4):27184. 48. Aarts MM, et al. The nucleolar targeting signal (NTS) of parathyroid hormone related protein mediates endocytosis and nucleolar translocation. J Bone Min Res. 1999;14(9):1493503. 49. Toribio RE, et al. The midregion, nuclear localization sequence, and C terminus of PTHrP regulate skeletal development, hematopoiesis, and survival in mice. FASEB J. 2010;24(6):194757. 50. Diefenbach-Jagger H, et al. Arg21 is the preferred kexin cleavage site in parathyroid-hormone-related protein. Eur J Biochem. 1995;229(1):918. 51. Imamura H, et al. Urinary excretion of parathyroid hormone-related protein fragments in patients with humoral hypercalcemia of malignancy and hypercalcemic tumor-bearing nude mice. J Bone Miner Res. 1991;6(1):7784. 52. Sirchia R, Luparello C. Mid-region parathyroid hormone-related protein (PTHrP) and gene expression of MDA-MB231 breast cancer cells. Biol Chem. 2007;388(5):45765. 53. Kumari R, Robertson JF, Watson SA. Nuclear targeting of a midregion PTHrP fragment is necessary for stimulating growth in breast cancer cells. Int J Cancer. 2006;119(1):4959. 54. Tran TH, et al. Loss of Nuclear localized parathyroid hormone-related protein in primary breast Cancer predicts poor clinical outcome and correlates with suppressed Stat5 signaling. Clin Cancer Res. 2018;24(24):635566. 55. Washam CL, et al. Identification of PTHrP(1248) as a plasma biomarker associated with breast cancer bone metastasis. Cancer Epidemiol Biomarkers Prev. 2013;22(5):97283. 56. Li X, et al. LIF promotes tumorigenesis and metastasis of breast cancer through the AKT-mTOR pathway. Oncotarget. 2014;5(3):788801. 57. Ghanei Z, et al. Immunization against leukemia inhibitory factor and its receptor suppresses tumor formation of breast cancer initiating cells in BALB/c mouse. Sci Rep. 2020;10(1):11465. 58. Viswanadhapalli S, et al. EC359: a first-in-class small-molecule inhibitor for Targeting Oncogenic LIFR Signaling in Triple-negative breast Cancer. Mol Cancer Ther. 2019;18(8):134154. 59. Bolin C, et al. Oncostatin m promotes mammary tumor metastasis to bone and osteolytic bone degradation. Genes Cancer. 2012;3(2):11730. 60. Wysoczynski M, et al. Leukemia inhibitory factor: a newly identified metastatic factor in Rhabdomyosarcomas. Cancer Res. 2007;67(5):213140. 61. Maruta S, et al. A role for leukemia inhibitory factor in melanoma-induced bone metastasis. Clin Exp Metastasis. 2008;26(2):133. 62. Ansari N et al. Autocrine and Paracrine Regulation of the murine skeleton by osteocyte-derived parathyroid hormone-related protein. 2018;33(1):13753. 63. Clements ME, Johnson RW. PREX1 drives spontaneous bone dissemination of ER + breast cancer cells. Oncogene. 2020;39(6):131834. 64. Johnson RW, et al. TGF-beta promotion of Gli2-induced expression of parathyroid hormone-related protein, an important osteolytic factor in bone metastasis, is independent of canonical hedgehog signaling. Cancer Res. 2011;71(3):82231. Publishers Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.  ...

... The item referenced in this repository content can be found by following the link on the descriptive page.  ...

... The item referenced in this repository content can be found by following the link on the descriptive page.  ...

... The item referenced in this repository content can be found by following the link on the descriptive page.  ...