Buscar

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

... Asian Spine Journal Asian Spine Journal Basic Study Asian Spine Spine J 2024;18(1):1-11 https://doi.org/10.31616/asj.2023.0203 Lumbar in Supine versus Weight-Bearing MRI 1 Lumbar Spine Anatomy in Supine versus WeightBearing Magnetic Resonance Imaging: Detecting Significant Positional Changes and Testing Reliability of Quantification Niladri Kumar Mahato1,2,3, Paramanand Maharaj4, Brian C. Clark2,3 1 College of Osteopathic Medicine, Marian University, Indianapolis, IN, USA 2 Ohio Musculoskeletal and Neurological Institute, Athens, OH, USA 3 Department of Biomedical Sciences, Ohio University, Athens, OH, USA 4 Department of Radiology, Eric Williams Medical Science Complex, University of the West Indies, St. Augustine, Trinidad and Tobago Study Design: Testing between and within group differences and assessing reliability of measurements. Purpose: To study and compare lumbar spine morphology in supine and weight-bearing (WB) magnetic resonance imaging (MRI). Overview of Literature: Upright lumbar MRI may uncover anatomical changes that may escape detection when using conventional supine imaging. This study quantified anatomical dimensions of the lumbar spine in the supine and WB MRI, compared specific morphometric differences between them, and tested the intra-rater reliability of the measurements. Repeated measures analysis was used to compare within- and between-session measurements performed on the supine and WB images. Reliability and agreement were assessed by calculating intraclass correlation (ICC) coefficient. Methods: Data from 12 adults without any history of back pain were used in this study. Sagittal T2-weighted images of the lumbar spine were acquired in the supine and WB positions twice (in two separate sessions scheduled within a week). Linear, angular dimensions, and cross-sectional areas (CSAs) were measured using proprietary software. Supine and WB data acquired from the two imaging sessions were tested for intra-rater reliability. Quantified data were normalized for each session to test the significance of differences. ICC was calculated to test the reliability of the measurements. Results: Linear, angular, and CSA measurements demonstrated strong within-position (supine and WB) correlations (r -values, 0.750.97). Between-position (supine vs. WB) differences were significant for all measured dimensions (p <0.05). Between-session measurements demonstrated a strong correlation (r -values, 0.640.83). Calculated ICC showed strong agreement among the measurements. Conclusions: Anatomical dimensions of the lumbar spine may demonstrate consistent and significant differences between supine and WB MRI for specific structural parameters. Keywords: Back pain; Sagittal; Spinal canal; Intervertebral; Degeneration Received Jun 30, 2023; Revised Aug 6, 2023; Accepted Aug 13, 2023 Corresponding author: Niladri Kumar Mahato College of Osteopathic Medicine, Evans Center, Marian University, Indianapolis, IN 46222, USA Tel: +1-812-603-3307, Fax: +1-317-955-6418, E-mail: nmahato@marian.edu ASJ Copyright 2024 by Korean Society of Spine Surgery This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited. Asian Spine Journal pISSN 1976-1902 eISSN 1976-7846 www.asianspinejournal.org 2 Niladri Kumar Mahato et al. Introduction Diagnostic magnetic resonance imaging (MRI) for low back pain (LBP) may not conclusively detect structural causes associated with the pain or dysfunction [1,2]. Studies have indicated that routine MRI, acquired in the supine position, could be inadequate or inconsistent in detecting structural changes causing mechanical or nonspecific back pain compared with imaging acquired in a weight-bearing (WB) position [3,4]. Conventional supine MRI used to diagnose LBP and spondylolisthetic segmental instability does not involve physiological weight bearing and may miss the detection of narrowing of the spinal canal, disks, and intervertebral foramina, as a cause of back pain, which could be identified upon imaging in WB positions [5,6]. Treatment decisions for disk herniation, canal stenosis, and intervertebral foramen narrowing are made based on clinical evaluation, which is supplemented by related radiological assessment. Accordingly, the effects of WB on diagnostic imaging may warrant further investigation [7,8]. Although WB MRI is not yet a standard for the diagnostic evaluation of LBP or listhetic segmental instability, instrumented spine loading in the supine and WB (sitting and standing) positions has been documented in experimental and diagnostic imaging studies [9-13]. Because MRI acquisition usually takes longer than X-ray techniques, MRI protocols could be physically challenging for patients with concurrent back pain. In contrast, spinal soft tissues are best visualized using MRI, and WB MRI has been demonstrated to be critical in detecting structural causes of LBP that may slip detection with imaging in the supine position [9-11]. Although some studies have examined the ability of positional MRI to uncover anomalous spine anatomy in patients with symptomatic lumbar stenosis, normative data comparing supine and WB lumbar spine anatomy in healthy individuals are limited, particularly in terms of back pain literature. Additionally, although positional MRI data comparing physiological effects of spine loading in healthy individuals are limited from single-session studies, tests for intra-rater reliability and statistical comparisons between positional data from multiple sessions are not readily available in the literature [8-12]. Our study used open WB MRI for lumbar spine imaging in a cohort of healthy participants (1) to detect the specific dimensions of the spine affected by positional changes in imaging, (2) to test whether WB induced sig- Asian Spine J 2024;18(1):1-11 nificant changes in these dimensions between the imaging positions, and (3) to evaluate the consistency of the within- and between-positional measurements for the detected dimensions. Materials and Methods This descriptive imaging study included a cohort of healthy participants to quantify morphometric changes in their lumbar spines induced by WB by comparing specific dimensions in supine versus WB MRI images. Additionally, we assessed the statistical significance of positional changes following WB and evaluated the consistency of our measurements using intra-rater reliability testing [13]. Imaging was performed in two separate sessions; each session involved imaging once in the WB position and once in the supine position. 1. Participant selection In this study, 13 volunteers aged 1860 years (mean: 3812 years; seven females and six males) participated. The exclusion criteria were as follows: individuals who reported LBP with a 1 on the numeric pain rating scale, those with a history of spinal surgery, those with any orthopedic or neurological impairment, those with cancer or tumors, those with cardiopulmonary disorders, those with clinical depression, and those who were taking medications or supplements for LBP. General compatibility for MRI was assessed; subjects with a body mass index (BMI) >32 kg/m2 or those with any physical limitations that impaired their ability to stand inside an open MRI system were excluded. Visit 1 included study orientation and subject consent. Visits 2 and 3 were scheduled 1 week apart at the university MRI facility [14]. One male participant dropped out of the study for personal reasons. Because of the relatively younger age of the study participants, none of the images showed any age-related degenerative changes that may have affected the anatomical parameters measured, except for two participants, who showed some evidence of age-related degeneration at their L3L4 and lumbosacral levels. Informed consent was obtained from all individual participants included in the study. The Institutional Review Board approved the study protocol (IRB #16-F-13), and all subjects gave written informed consent. Asian Spine Journal Lumbar Spine in Supine versus Weight-Bearing MRI 3 2. Imaging protocol Scout images were obtained iso-centered at the L3L4 intervertebral disks after the participants were secured in the supine position on the MRI table. The height of the footrest was adjusted and recorded in the supine position before the table was slowly rotated to 87 vertical (as specified by the machine-operating protocol to prevent fall) to the ground, as the participants eased into a WB position with their feet shoulder-width apart. The participants were comfortably secured above the knee and at the pelvis using cushioned straps to help them maintain an upright position during imaging. All images were acquired using an open-field magnet (0.25 T) G Scan Brio WB MRI (Esaote S.p.A., Genova, Italy) in an upright position. Next, the participants were placed in the supine position by tilting the table back to the horizontal plane. They laid supine and rested for 5 minutes to relax their spine before the initiation of image acquisition in the supine position. Sagittal MR images of the lumbar spine (L1S1) were acquired in both positions using a fast spin-echo T2 sequence (TR=3,520 ms; TE=125 ms; number of acquisitions=1; matrix=288234; FOV=320320; oversampling=185%; slice thickness=4 mm; gap=1 mm; acquisition time=4 minutes 41 seconds). The scanning time and sequences we followed have been described in earlier studies reporting <5 minutes of scan time at each position [15]. Typically, each participant spent approximately 45 minutes per session, which included positioning, tilting, image acquisition, and resting time. The height of the standing platform, positions of the cushion straps, and the lumbar coil placements were measured for each participant to be replicated in the subsequent imaging sessions. 3. Image analysis All Digital Imaging and Communications in Medicine images were transferred to a proprietary image analysis platform (OrthoCAD; Esaote S.p.A.). The borders of the lumbar spine, first sacral vertebral bodies, spinal canal, and vertebral foramina were manually outlined on the image slices using the software. Images were then segmented and digitized using proprietary three-dimensional (3D) segmentation software (OrthoCAD) to quantify anatomical dimensions. The first author (N.K.M.) used semiautomated segmentation and measurement tools available in the software to quantify all dimensions presented in this Fig. 1. Mid-sagittal lumbar spine T2-weighted image showing the sacral angle, disc heights, and mid zonal spinal canal dimensions as demarcated by the proprietary software tool OrthoCAD. study (Fig. 1). Additionally, a single observer used ImageJ (National Institutes of Health, Bethesda, MD, USA; Imagej.nih.gov/ij/docs/faqs.html) to verify the linear dimensions measured in the images. 4. Measurement variables Linear dimensions (in millimeters) were measured as follows: Sagittal dimensions of the spinal canal were measured as the minimal available space at the vertebral junctional areas, referred to as zones (Fig. 1). Disk heights were measured in the mid-intervertebral spaces. Foraminal heights were measured at the narrowest portions along the intervertebral openings. All segmental and intersegmental angles (in degrees) were measured in the sagittal plane using the endplates superior to the concerned vertebrae. The foramen cross-sectional area (CSA) (in square millimeters) at each level was measured at the narrowest available point along image slices capturing the intervertebral foramen. Slice positions and imaging parameters in the protocol were kept consistent between the two sessions. 5. Statistical analysis Supine and WB measurements obtained from the two 4 Niladri Kumar Mahato et al. Asian Spine J 2024;18(1):1-11 sessions were tested separately for correlation using Pearsons correlation coefficient. Within-session differences between supine and WB images were tested for statistical significance using a two-tailed paired Student t-test. Considering potential changes in magnification and minor variability associated with participant positioning within the scanner in the two imaging sessions, the differences calculated between the supine and WB images were normalized for each session. This was performed by subtracting the supine measurements from the WB measurements and then dividing that value by the WB measurement (at each level) to obtain an absolute value for the dimension: [(Supine-WB)/WB]. Correlations between normalized differences between the two sessions were determined using Pearsons correlation coefficient. Additionally, factorial repeated measured analysis of variance (ANOVA) (without replication; 2 positions2 sessions), followed by multiple post-hoc comparisons (Tukeys test), was performed to determine the statistical significance of differences in each dimension across images between both positions and both imaging sessions. The intra-rater reliability of the measurements was evaluated by computing the intraclass correlation coefficients (ICC) of all measurement variables for the two sessions and two scanning positions (IBM SPSS Statistics for Windows, ver. 24.0; IBM Corp., Armonk, NY, USA). Results 1. Linear dimensions 1) Spinal canal dimensions Zones increased with WB (Fig. 2). Our findings showed that this increase was proportional to the initial dimensions measured across each zone in the supine position. We observed a strong correlation and intra-rater reliability for measurements acquired in the supine and WB positions (r=0.89 and 0.96, respectively) (Table 1). Furthermore, between-session normalized differences at the five zones were strongly correlated, indicating consistency in detecting supine versus WB changes (r=0.83). Additionally, differences in sagittal dimensions between the supine and WB positions were statistically significant in both sessions (F [3,59]=6.61; p<0.01). The spinal canal dimensions showed an increase at all levels from the supine position to the WB position. The change percentages observed were as follows: L1L2=11.35%, L2L3=10.18%, L3L4=3.21%, L4L5=8.93%, and L5S1=12.56%. 2) Disk heights Disk heights (mm) were reduced with WB (Fig. 2). Decreases in this dimension corresponded proportionally to the initial heights at each level in the supine position [15]. We detected a strong correlation and intra-rater reliability in the measurements acquired in the supine and WB positions (r=0.94 and 0.91, respectively) (Table 1). A moderate correlation was detected in the between-session normalized differences in disk heights at all five intervertebral Spinal canal thickness and disc height Supine Weight bearing 20 18 16 14 12 10 8 6 4 2 0 L1L2 L2L3 L3L4 L4L5 Spinal canal thickness L5S1 L1L2 L2L3 L3L4 L4L5 L5S1 Disc height Fig. 2. Spinal canal dimensions and disc heights (in millimeters) measured by supine and weight-bearing imaging. All parameters returned significant values when comparing the two imaging positions (p <0.05). Asian Spine Journal Lumbar Spine in Supine versus Weight-Bearing MRI 5 Table 1. Descriptive statistics of linear dimensions measured in images acquired with supine and WB positions in two separate sessions Zone Variable L1L2 L2L3 L3L4 L4L5 L5S1 r-value Spinal canal sagittal dimension Supine sessions (mm) 13.951.34 13.391.38 13.501.15 12.861.23 12.612.00 0.89 WB sessions (mm) 15.541.22 14.761.46 13.941.35 14.011.31 14.202.61 0.96 11.35 10.18 3.21 8.93 12.56 0.83 Supine sessions (mm) 11.172.13 11.631.30 13.441.74 14.433.10 12.971.71 0.94 WB sessions (mm) 10.111.42 11.261.57 12.841.62 13.812.66 12.142.31 0.91 9.50 3.20 4.44 4.41 6.40 0.64 Supine sessions (mm) 16.933.45 17.974.07 17.873.58 16.633.78 14.322.32 0.96 WB sessions (mm) 15.932.77 16.923.81 17.483.06 15.953.35 15.392.71 0.97 12.15 7.63 7.87 4.07 7.49 0.74 Supine sessions (mm) 16.663.48 17.194.22 18.865.05 13.894.29 14.493.41 0.96 WB sessions (mm) 14.632.36 15.883.11 17.383.51 15.682.66 14.72.21 0.97 5.88 5.86 2.18 5.93 2.47 0.79 a) % Changes Disc height a) % Changes Intervertebral foramen height Left % Changesa) Right a) % Changes Values are presented as meanstandard deviation or %. r-values indicate Pearsons correlation coefficients. Percentage changes denote measurement changes from the supine to the weight bearing positions. Percentage changes with indicate a decrease. WB, weight bearing. a) Supine vs. WB differences (normalized): normalized supine vs. WB difference for each session=(supineWB)/WB. levels, indicating consistency in detecting supine versus WB changes (r=0.64). Differences in disk heights between the supine and WB positions were statistically significant for both sessions (F [3,59]=23.34; p<0.01). The disk heights decreased at all levels from the supine position to the WB position. The change percentages observed were as follows: L1L2=9.50%, L2L3=3.20%, L3L4=4.44%, L4L5=4.41%, and L5S1=6.40%. 3) Foraminal heights Foraminal heights (mm) were bilaterally reduced with WB, except at the L5S1 junctions, where the foramina appeared to slightly enlarge. Analysis revealed a strong correlation and intra-rater reliability in the measurements acquired in the supine and WB positions bilaterally (r=0.96 and 0.97, respectively) (Table 1). The betweensession normalized differences in foraminal heights were strongly correlated at all five intervertebral levels, indicating consistency in detecting supine versus WB changes (right: r=0.79; left: r=0.74). The differences in foraminal heights between the supine and WB positions were sta- tistically significant in both sessions (F [3,59]=20.75; p<0.01). The intra-rater reliability of the measurements was 75% for each session and each scan position. All foraminal heights showed an increase at all levels, on both sides, from the supine position to the WB position, except for L5S1, which showed a decrease bilaterally. The change percentages observed on both sides were as follows: right: L1L2=12.15%, L2L3=7.63%, L3L4=7.87%, L4L5=4.07%, and L5S1=7.49%; left: L1L2=5.88%, L2L3=5.86%, L3L4=2.18%, L4L5=5.93%, and L5 S1=1.47%. 2. Angular dimensions: segmental angles Segmental angles increased with WB. In our findings, the angles increased caudally and were proportionate to the magnitude of the baseline angles measured in the supine position. However, the lumbosacral angle was reduced with WB (21.964.55 versus 20.794.99). Among the intersegmental angles, the lordosis (L1S1), sacral (horizontalS1), lumbar (L1L5), and lower lordosis (L3 6 Niladri Kumar Mahato et al. Asian Spine J 2024;18(1):1-11 70 Angular dimension Supine Weight bearing 60 50 40 30 20 10 0 L1L2 L2L3 L3L4 L4L5 L5S1 L3L4 angle L1L5 angle Sacral angle L1S1 angle Sehmental Inter-segmental Fig. 3. Segmental and inter-segmental angles (in degrees) measured by supine and weight-bearing imaging. All parameters returned significant values when comparing the two imaging positions (p <0.05). Table 2. Descriptive statistics of angular dimensions measured in images acquired in supine and WB positions in two separate sessions Segmental angles L1L2 Supine sessions () 3.331.24 WB sessions () 5.251.54 57.50 41.60 a) % Changes L2L3 L3L4 L4L5 L5S1 r -value 5.212.23 7.543.35 7.382.16 11.632.26 12.925.69 21.964.55 0.91 15.544.47 20.794.99 0.97 54.14 20.32 6.45 0.74 Values are presented as meanstandard deviation or %. r -values indicate Pearsons correlation coefficients. Percentage changes denote measurement changes from the supine to the weight bearing positions. Percentage changes with indicate a decrease. WB, weight bearing. a) Between-sessions differences (normalized): normalized supine vs. WB difference for each session=(supineWB)/WB. Table 3. Descriptive statistics of angular dimensions measured in images acquired in supine and WB positions in two separate sessions Inter-segmental angle Lordosis angle L1S1 Lumbar angle L1L5 Sacral angle horizontalS1 Lower lordosis angle L3L5 MeanSD r -value MeanSD r -value MeanSD r -value MeanSD r -value Supine sessions () 47.679.24 0.98 27.049.39 0.94 36.387.20 0.75 22.006.63 0.93 WB sessions () 56.337.67 0.94 38.585.42 0.88 37.678.51 0.78 27.255.46 0.76 23.86 0.84 39.29 0.54 6.07 0.61 18.18 0.68 a) % Changes r -values indicate Pearsons correlation coefficients. Percentage changes denote measurement changes from the supine to the weight bearing positions. Percentage changes with indicate a decrease. WB, weight bearing; SD, standard deviation. a) Supine vs. WB differences (normalized): normalized supine vs. WB difference for each session=(supineWB)/WB. L5) angles were the largest to the smallest in the supine position (Fig. 3). These angles increased with WB, with the highest increase observed in the lumbar angle (L1 L5) and the lowest increase observed in the sacral angle (horizontalS1) while standing. The measurements acquired in the supine and WB positions exhibited a strong correlation and intra-rater reliability (r=0.91 and 0.97, respectively) (Tables 2, 3). A moderately strong correlation was detected for normalized differences in between- session segmental angular data, indicating consistency in detecting supine versus WB alterations (r=0.74). Angular parameters were significantly different between the supine and WB positions in both sessions (F [3,59]=72.16; p<0.01). The intra-rater reliability of the measurements was 75% for each session and each scan position. All angular dimensions decreased at all levels from the supine position to the WB position, except for the L5 S1 intervertebral segmental angle, which increased. The Asian Spine Journal Lumbar Spine in Supine versus Weight-Bearing MRI 7 Foramen CSA Supine Weight bearing 160 140 120 100 80 60 40 20 0 L1L2 L2L3 L3L4 L4L5 L5S1 L1L2 L2L3 L3L4 L4L5 L5S1 Left Right Fig. 4. Intervertebral foramen cross-sectional areas (CSAs) measured by supine and weight-bearing imaging. All parameters returned significant values when comparing the two imaging positions (p <0.05). Table 4. Descriptive statistics of CSAs measured in images acquired in the supine and WB positions in two separate sessions L1L2 L2L3 L3L4 L4L5 L5S1 r -value 118.8719.79 114.0323.68 113.4915.54 116.1818.43 124.2516.16 0.87 109.2112.26 108.6019.25 108.5215.58 109.3316.60 125.2013.22 0.97 8.13 4.76 4.38 5.90 1.28 0.68 117.7015.64 112.4917.98 109.0315.18 110.73147.48 114.35129.32 0.92 110.7915.88 110.1312.21 111.3914.29 109.3615.88 122.2613.22 0.97 5.86 2.10 2.16 1.24 6.92 0.74 Inter-vertebral foramen CSA Left Supine sessions (mm2) 2 WB sessions (mm ) a) % Changes Right Supine sessions (mm2) 2 WB sessions (mm ) % Changesa) Values are presented as meanstandard deviation or %. r -values indicate Pearsons correlation coefficients. Percentage changes denote measurement changes from the supine to the weight bearing positions. Percentage changes with indicate a decrease. CSA, cross-sectional area; WB, weight bearing. a) Supine vs. WB differences (normalized): normalization of supine vs. WB difference for each session=(supineWB)/WB. change percentages observed were as follows: intervertebral segmental angles: L1L2=57.50%, L2L3=41.60%, L3L4=54.14%, L4L5=20.32%, L5S1=6.45%; L3L5 angle=23.86%; L1L5 angle=39.29%; sacral angle horizontalS1=6.07%; and L1S1 angle=18.18%. 3. Intervertebral foraminal CSA Foraminal CSAs (mm2) decreased with WB, except at the lumbosacral junction (L5S1), which consistently increased bilaterally (left: 124.8516.16 versus 125.2013.22; right: 114.35129.32 versus 122.2613.22, respectively) (Fig. 4). Increased CSA was also noted at the mid-lumbar level (L3L4) on the right side. A strong correlation and intra-rater reliability were observed in the CSAs acquired in the supine and WB positions (left: r=0.87 and 0.97; right: r=0.92 and 0.97, respectively) (Table 4). Furthermore, a strong correlation between the normalized differences from the two sessions was detected at all five levels, indicating consistency in measuring changes induced by WB (left: r=0.68; right: r=0.74). Foraminal CSAs significantly differed between the supine and WB positions in both sessions (F [3,59]=5.59; p<0.01). The intra-rater reliability of the measurements was 75% for each session and each scan position. The ICC showed high values of agreement for all dimensions (measured by sessions and scanning positions) calculated as inter-rater reliability using absolute agreement and a two-way random model 8 Niladri Kumar Mahato et al. Asian Spine J 2024;18(1):1-11 Table 5. ICC calculated for measured variables by sessions and by scanning positions Supine Variable WB ICC 95% CI ICC 95% CI Sagittal dimension 0.72 0.580.71 0.77 0.650.85 Disc heights 0.88 0.810.93 0.91 0.860.94 Left 0.94 0.910.96 0.94 0.910.96 Right 0.95 0.910.97 0.96 0.930.96 0.97 0.960.98 0.98 0.970.98 Left 0.92 0.870.95 0.83 0.740.90 Right 0.84 0.750.90 0.89 0.820.93 Linear dimensions Intervertebral foramen height Angular dimensions Intervertebral foramen CSA Inter-rater reliability was performed using absolute agreement and a two-way random model reporting single measures for ICC, in SPSS. ICC, intraclass correlation coefficients; WB, weight bearing; CI, confidence interval; CSA, cross-sectional area. (Table 5). All intervertebral foraminal CSAs increased at all levels, on both sides, from the supine position to the WB position, except for the L5S1 level, which decreased bilaterally. The change percentages observed on both sides were as follows: right: L1L2=8.13%, L2L3=4.76%, L3L4=4.38%, L4L5=5.90%, L5S1=0.28%; left: L1 L2=5.86%, L2L3=2.10%, L3L4=2.16%, L4L5=1.24%, and L5S1=6.92%. Only one pairwise Tukey honestly significant difference comparison (i.e., foraminal height) between the supine and WB positions was significant (p<0.01); all other comparisons were insignificant. This result is expected considering that the pairs of the compared segmental means (two supine and two WB) remained clustered together. All positional comparisons tested (and presented) returned significant p-values (p<0.05) when assessed using paired Student t -test; this result was supported by the significant ANOVA F statistic, as reported above. Discussion This study assessed (1) specific dimensions that demonstrated morphometric shifts between the two positions of imaging, (2) the reliability of measurements on lumbar spine images between two loading positions and across two imaging sessions, and (3) the statistical significance of the anatomical changes observed between the supine and WB positions. Increases in the sagittal dimensions observed at the vertebral junctional zones within the spinal canal may have occurred because of dynamic alterations in intra- and intersegmental changes in the vertebral angles, increased apposition of the zygapophyseal joints, and bulging of the interlaminar ligament system with WB, as reported in the literature [16,17]. Such dynamic effects on the lumbar spinal canal with axially loaded spine observed with computed tomographymyelography and MRI in patients with sciatica have been reported and agree to the findings of this study involving asymptomatic individuals [6,11]. These increments could also reflect increasing curvature of the lumbar spine, particularly at the thoracolumbar and lumbosacral junctions with the upright spine [18-21]. Reductions in disk heights reflect axial loading, though the reductions are not proportionately as large as reported in some earlier studies [22-24]. In this study, disk heights were measured at the mid-intervertebral spaces to detect changes at the level of the centrally placed nucleus. With upright WB, compression of the peripheral annulus may stabilize or even increase the height of the centrally located disk nucleus. Disk heights were reduced at the L1 L2 junction. Additionally, although earlier studies have reported Modic I changes to increase in size from the supine position to the standing position with MRI studies that correlated with an increase in pain intensity in the standing position, this study further substantiates the evidence that WB may uncover the biomechanical stress and active discopathy theories in LBP, as suggested by the authors [15]. Intervertebral foraminal heights decreased upon standing, possibly due to compressive loading at the posterior vertebral elements from the accentuation of the overall lumbar curvature, reduction of peripheral disk heights, or increased anterior angulation at segmental endplates [19]. This aligns with the observation that forward flexion facilitates foraminal widening by reducing posterior element stress, increasing the vertical dimension of the foramen [22-24]. Accordingly, forward flexion of the lumbar spine often mechanically reduces pain resulting from foraminal and spinal canal narrowing [25-27]. Interestingly, the L5 S1 foramina marginally increased in this dimension on WB, possibly due to a greater increase in the sacral angle upon standing. The relatively coronal orientation of the facet joints at this level may serve as struts that not only prevented the L5 listhesis at this junction but also acted as Asian Spine Journal a fulcrum to help increase the L5S1 angle, thereby opening up the foramen dimensions [20-22]. Interestingly, some studies have reported comparable changes in the intervertebral foraminal spaces as effects of axially loaded MRI of the lumbar spine on dural sac and lateral recesses [16]. All segmental and intersegmental angles in the sagittal plane increased with different magnitudes with WB. However, the overall L5S1 angle slightly decreased (with greater variance). Nonetheless, all other angles related to the S1 endplate, including the lordosis and sacral angles, significantly increased with the upright position [23,24]. Further investigation is required to detect whether lumbar spine curvature variability is associated with overall changes in the spinal canal or intervertebral foraminal dimensions in individuals with back pain upon WB [25]. Interestingly, foraminal CSAs demonstrated changes due to alterations in disk and foraminal heights on WB. However, the lumbosacral foramina marginally increased bilaterally in CSA on WB. Thus, unlike the rest of the cranially located foramina, the L5S1 foramina widened on standing, with a concomitant increase in CSA and height. Intuitively, one may assume that although the vertical foraminal heights reduce on WB, the available anteroposterior space in these foramina may increase, thereby preventing critical narrowing of the passage [26]. However, this study indicates that the corresponding CSAs at the L1L4 foramina also decrease with standing. Notwithstanding, the lumbar intervertebral foramina are tubular spaces that are directed somewhat posterolaterally (not strictly in the coronal plane). Therefore, any change in the dimension perpendicular to the plane of the sagittal image may remain undetected and could prevent the detection of an increase in the foraminal dimension in that plane [27]. The CSA at the lumbosacral (L5S1) foramen marginally increased bilaterally, as did its vertical dimensions. This increase may be due to the significant angular increases associated with the S1 endplate on WB. Because narrowing of the foramina is pathologically related to nerve compression and sciatica, changes in all foraminal CSAs on WB could yield clinically valuable information [25,28]. This study has some limitations. First, none of the study participants had chronic or acute back pain episodes, and all individuals were within a specific BMI range. Thus, we assume that the observed anatomical changes were within physiological limits and represent the ef- Lumbar Spine in Supine versus Weight-Bearing MRI 9 fects of normal WB in the upright position. Second, the disk hydration status was not objectively accounted for every imaging session in terms of its potential effects on the disk height measurements. The disk height changes reported in this study were discussed as changes that incurred with gravitational load bearing only. Some changes related to disk hydration may have impacted the variability of data reported in this study. Lastly, although the scan time in the upright position was kept to the minimal mainly to acquire sagittal images and very few axial images for comparison, shortening the scan time further in the protocol may help include participants with back pain in such studies in the future. Comparisons could also be improved by generating a 3D rendition of the images and measuring additional parameters in multiaxial planes. Examining specific anatomical dimensions that change in the WB position could be critical in identifying specific structural causes of LBP and lumbar instability, for example, markers that may remain undetected on supine MRI [23]. Despite the relatively lower sample size, the number of images examined could delineate (normality of distributions, statistical significance within the assigned alpha levels) differences in the anatomical parameters measured in the study. Images acquired from the relatively low-strength 0.25T MRI machine could provide adequate resolution to allow accurate and consistent measurements. The equipment is approved for clinical usage and comes with validated measurement software for use in clinical research and practice. Measurements of the anatomical changes presented in this study may help contextualize the importance of WB MRI in detecting potential anatomical changes that may be etiologically associated with back pain and spine instability. These changes may otherwise remain undetected with routine supine imaging. Though longer MRI protocols may be challenging for individuals with concurrent back pain, WB imaging could still be useful in detecting diagnostically valuable anatomical changes in the lumbar spine where indicated [29,30]. Conclusions WB imaging may improve detection of changes in the lumbar spine anatomy compared with supine MRI. As a future step, findings from this study could be further corroborated in individuals with back pain to assess the validity of assessment in the context of specific dimen- 10 Niladri Kumar Mahato et al. sions (concerning spinal canal, disks, and intervertebral foramina) that may be pathologically related to back pain. Conflict of Interest No potential conflict of interest relevant to this article was reported. ORCID Niladri Kumar Mahato: https://orcid.org/0000-0001-54391172; Paramanand Maharaj: https://orcid.org/0000-00031859-1606; Brian C. Clark: https://orcid.org/0000-00026021-6431 Author Contributions Conceptualization: NKM, BCC; data curation: NKM, BCC; formal analysis: NKM, BCC; funding acquisition: NA; methodology: NKM, BCC, PM; project administration: BCC, NKM; visualization: NKM; writingoriginal draft: NKM, BCC; writingreview & editing: NKM, BCC, PM; final approval of the manuscript: all authors. References 1. Weishaupt D, Schmid MR, Zanetti M, et al. Positional MR imaging of the lumbar spine: does it demonstrate nerve root compromise not visible at conventional MR imaging? Radiology 2000;215:247-53. 2. Baker MA, MacKay S. Please be upstanding: a narrative review of evidence comparing upright to supine lumbar spine MRI. Radiography (Lond) 2021;27:721-6. 3. Ozawa H, Kanno H, Koizumi Y, et al. Dynamic changes in the dural sac cross-sectional area on axial loaded MR imaging: is there a difference between degenerative spondylolisthesis and spinal stenosis? AJNR Am J Neuroradiol 2012;33:1191-7. 4. Maher C, Underwood M, Buchbinder R. Non-specific low back pain. Lancet 2017;389:736-47. 5. Mahato NK, Sybert D, Law T, Clark B. Effects of spine loading in a patient with post-decompression lumbar disc herniation: observations using an open weightbearing MRI. Eur Spine J 2017;26(Suppl 1):17-23. 6. Tarantino U, Fanucci E, Iundusi R, et al. Lumbar spine MRI in upright position for diagnosing acute and chronic low back pain: statistical analysis of morpho- Asian Spine J 2024;18(1):1-11 logical changes. J Orthop Traumatol 2013;14:15-22. 7. Jinkins JR, Dworkin JS, Green CA, et al. Upright, weight-bearing, dynamic-kinetic MRI of the spine pMRI/kMRI. Riv Neuroradiol 2002;15:333-57. 8. Kanno H, Endo T, Ozawa H, et al. Axial loading during magnetic resonance imaging in patients with lumbar spinal canal stenosis: does it reproduce the positional change of the dural sac detected by upright myelography? Spine (Phila Pa 1976) 2012;37:E985-92. 9. Wang J, Yang X. Age-related changes in the orientation of lumbar facet joints. Spine (Phila Pa 1976) 2009;34:E596-8. 10. Wessberg P, Danielson BI, Willen J. Comparison of Cobb angles in idiopathic scoliosis on standing radiographs and supine axially loaded MRI. Spine (Phila Pa 1976) 2006;31:3039-44. 11. Willen J, Danielson B, Gaulitz A, Niklason T, Schonstrom N, Hansson T. Dynamic effects on the lumbar spinal canal: axially loaded CT-myelography and MRI in patients with sciatica and/or neurogenic claudication. Spine (Phila Pa 1976) 1997;22:2968-76. 12. Kanno H, Ozawa H, Koizumi Y, et al. Changes in lumbar spondylolisthesis on axial-loaded MRI: do they reproduce the positional changes in the degree of olisthesis observed on X-ray images in the standing position? Spine J 2015;15:1255-62. 13. Kottner J, Audige L, Brorson S, et al. Guidelines for Reporting Reliability and Agreement Studies (GRRAS) were proposed. J Clin Epidemiol 2011;64:96-106. 14. Mahato NK, Montuelle S, Clark BC. Assessment of in vivo lumbar inter-vertebral motion: reliability of a novel dynamic weight-bearing magnetic resonance imaging technique using a side-bending task. Asian Spine J 2019;13:377-85. 15. Splendiani A, Bruno F, Marsecano C, et al. Modic I changes size increase from supine to standing MRI correlates with increase in pain intensity in standing position: uncovering the biomechanical stress and active discopathy theories in low back pain. Eur Spine J 2019;28:983-92. 16. Bulja D, Strika J, Jusufbegovic M, et al. Effects of axial loaded magnetic resonance imaging of lumbar spine on dural sac and lateral recesses. J Health Sci 2021;11:181-5. 17. McKay G, Torrie PA, Bertram W, et al. Myelography in the assessment of degenerative lumbar scoliosis Asian Spine Journal 18. 19. 20. 21. 22. 23. 24. and its influence on surgical management. Korean J Spine 2017;14:133-8. Jayakumar P, Nnadi C, Saifuddin A, Macsweeney E, Casey A. Dynamic degenerative lumbar spondylolisthesis: diagnosis with axial loaded magnetic resonance imaging. Spine (Phila Pa 1976) 2006;31:E298-301. Nguyen HS, Doan N, Shabani S, et al. Upright magnetic resonance imaging of the lumbar spine: back pain and radiculopathy. J Craniovertebr Junction Spine 2016;7:31-7. Chan DD, Gossett PC, Butz KD, Nauman EA, Neu CP. Comparison of intervertebral disc displacements measured under applied loading with MRI at 3.0 T and 9.4 T. J Biomech 2014;47:2801-6. Liu X, Zhao X, Long Y, et al. Facet sagittal orientation: possible role in the pathology of degenerative lumbar spinal stenosis. Spine (Phila Pa 1976) 2018;43:955-8. Mahato NK. Association of rudimentary sacral zygapophyseal facets and accessory and ligamentous articulations: implications for load transmission at the L5-S1 junction. Clin Anat 2010;23:707-11. Meakin JR, Gregory JS, Smith FW, Gilbert FJ, Aspden RM. Characterizing the shape of the lumbar spine using an active shape model: reliability and precision of the method. Spine (Phila Pa 1976) 2008;33:807-13. Meakin JR, Smith FW, Gilbert FJ, Aspden RM. The effect of axial load on the sagittal plane curvature of the upright human spine in vivo. J Biomech 2008;41:2850-4. Lumbar Spine in Supine versus Weight-Bearing MRI 11 25. Tebet MA. Current concepts on the sagittal balance and classification of spondylolysis and spondylolisthesis. Rev Bras Ortop 2014;49:3-12. 26. Meakin JR, Gregory JS, Aspden RM, Smith FW, Gilbert FJ. The intrinsic shape of the human lumbar spine in the supine, standing and sitting postures: characterization using an active shape model. J Anat 2009;215:206-11. 27. Keorochana G, Taghavi CE, Lee KB, et al. Effect of sagittal alignment on kinematic changes and degree of disc degeneration in the lumbar spine: an analysis using positional MRI. Spine (Phila Pa 1976) 2011;36:893-8. 28. Hebelka H, Rydberg N, Hutchins J, Lagerstrand K, Brisby H. Axial loading during MRI induces lumbar foraminal area changes and has the potential to improve diagnostics of nerve root compromise. J Clin Med 2022;11:2122. 29. Nordberg CL, Hansen BB, Nybing JD, et al. Weightbearing MRI of the lumbar spine: technical aspects. Semin Musculoskelet Radiol 2019;23:609-20. 30. Charoensuk J, Laothamatas J, Sungkarat W, Worapruekjaru L, Hooncharoen B, Chousangsuntorn K. Axial loading during supine MRI for improved assessment of lumbar spine: comparison with standing MRI. Acta Radiol 2023;64:217-27.  ...

... 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  ...

... Indiana Audubon Society NEWSLETTER FEB.MARCH 2024 Vol.55, No.1 CLAY-COLORED SPARROW DISCOVERY Indiana's first confirmed breeding record A LOWER 48 BIG YEAR A Hoosier professor's 2023 sabbatical adventure MEMBERSHIP PERKS Electronic bi-monthly & quarterly newsletters Mailed copies of bi-monthly & quarterly newsletters Access to Cornell Lab's Birds of the World ($49/yr. value) Early registration to Indiana Dunes Birding Festival Discounted Field Trips & Programs Special Access to Mary Gray Birding Sanctuary 15% discount on all IAS online store merchandise Register a friend to any IAS field trip at member rate 0 UPCOMING Join Indiana Audubon for an assortment of field trips and workshops being offered this summer and early fall. Visit indianaaudubon.org/events for more information or to register for any of these upcoming activities. Additional events are listed on our website, and many fill early! Register online today. FEB. 13: Feathered Focus: Bird Safe Indy [Virtual, Zoom] FEB. 16: Tropical Birding Live from Amagusa Reserve, Ecuador [Virtual, Zoom] FEB. 23: Hoosier Birders' Hour: Feeding Backyard Birds [Virtual, Zoom] Clay-colored Sparrow by Ryan Sanderson MARCH 12: Feathered Focus: The Secret to Attracting Birds [Virtual, Zoom] MARCH 16: Flock & Learn: Waterfowl Walk at Eagle Creek Park [Indianapolis] MARCH 16: Summit Lake Spring Waterfowl Trip [New Castle] MARCH 16: "What's That Waterfowl?" Young Birders Walk [Indianapolis] MARCH 23: Kankakee Ducks and More Early Spring Birding [North Judson] MARCH 26: Indiana Birding: Crash Course [Virtual, Zoom] 1 CARDINAL | INDIANA AUDUBON SOCIETY BN EG LL Be a part of Indiana Audubon! Sign up at: Indianaaudubon.org/membership TRIPS & EVENTS On the cover: E S TU S RK S 3 -1 UD E ST AG PA A BE T/I Y EN HI RS NE RA RT LI B CL CO EX 8 O R C O RY P W I TH PE $1 0 R W ITH IFE ME 0 IA S BE & ES IV E$ 7 A 75 L ILE AL E P R IN IV DE $2 NT ND YO UB S C G L R B E RD FA TR Y B E JO 0 $1 0 0 $5 FIT G NE K C T TI N O OC N E IBU 0 LY FL MI $4 E F A CO EN N TH ILY N UP FA M SIG US 0 LT $3 R 1 AD U FO I DU FIT AL S NE BE DIV TODAY LI F The Cardinal Newsletter is a bi-monthly publication of the Indiana Audubon Society. Its purpose is to share stories and conversations so that members and the birding community beyond can stay meaningfully connected both to birds and to the people dedicated to their protection. RENEW CA Contributors Leah Baker David Benson Wesley Homoya Joni James Libby Keyes Theresa Murray Ryan Sanderson Kristin Stratton OR SIC Editor & Graphic Designer Whitney Yoerger JOIN IN Executive Director Brad Bumgardner BA Production Team NEWS INDIANA AUDUBON SOCIE T Y W E'R E H I R I N G! Indiana Audubon is in search of a dedicated individual to become part of our staff as the Development and Engagement Manager. This role is crucial for advancing Indiana Audubons mission and involves actively contributing to fundraising and development initiatives, ensuring the ongoing success of our conservation, education, and research programs and initiatives. The job posting will be open until Feb. 15 or until the position is filled. Apply today! Get the details at indianaaudubon. org/currentopenings or scan the QR code. Trout-lily at Mary Gray Bird Sanctuary by Libby Keyes. The Sanctuary entrance by Libby Keyes. H E A LTH Y H I K I N G I N 2024 This year, immerse yourself in nature and improved well-being by exploring the 8 miles of trails at Mary Gray Bird Sanctuary in Connersville. Delight in the emotional and physical benefits of the sanctuary as you aim to cover a minimum of 25 miles there this year. Relish the peace, encounter birds and other wildlife, and experience this Fayette County gem while walking, hiking, or running on the trails. Interested in joining the MGBS Healthy Hiking Challenge? The sanctuary welcomes visitors from sunrise to sunset every day of the year. And feel free to bring your family and friends along for the adventure. Trail maps are available at the kiosk near the shelter across from the centralized parking lot. While there is no entrance fee to visit, donations are welcome and will help keep the sanctuary open to all visitors throughout the years to come. The sanctuary is not a state park, receives no federal or state funding. Visitor and member donations help keep the Sanctuary the beautiful place that it is. Take your first step by downloading the MGBS Healthy Hiking Challenge Entry Form at indianaaudubon.org/hike. FEB. MARCH 2024 2 NEWS INDIANA AUDUBON SOCIE T Y PR E PA R E FO R TH E D U N E S B I R D I N G FE S TI VA L Get ready for an even bigger and better Indiana Dunes Birding Festival, taking place from May 16-19! Featuring an exciting lineup of over 180 field trips, programs, and workshops, our festival schedule promises an exhilarating birding experience. Explore the full schedule and plan your festival activities in advance. Easily peruse the schedule or access the guidebook PDF online. New this year, Indiana Audubon members must renew their membership by midnight on Feb. 28 to secure early access to festival registration. Priority registration for Indiana Audubon members begins at 10 a.m. Central Time on March 1. Secure your spot for the four-day event at just $65, with options available for youth and single-day participation. We look forward to celebrating the festival's 10th anniversary with you! Make your festival plans, mark your calendar for priority registration on March 1, and renew your membership by midnight on Feb. 28. View the schedule at indunesbirdingfestival.com/pages/schedule or scan the QR code on the right. B R I DG I N G V I S I O N A N D AC TI O N Indiana Audubon's annual winter retreat, held from January 12 to 14 in Rockville, proved to be a productive weekend for the IAS Board of Directors. Despite the freezing temperatures, the board engaged in fruitful planning sessions, outlining the path for numerous exciting programs, committee tasks, and initiatives at Mary Gray Bird Sanctuary and for the organization as a whole, ensuring a promising and eventful 2024. Board members also took moments to enjoy birding from the windows. A sincere THANK YOU goes to our dedicated Board of Directors and staff, whose commitment made this weekend especially impactful. 3 CARDINAL | INDIANA AUDUBON SOCIETY Grant Report B E THAN Y L AN D I N STITUTE BY THERESA MURR AY On April 9, 2022, the Avian Conservation Uganda Society (ACUS) visited Bethany Land Institute (BLI) in Luwero District, Uganda, as part of a pilot program. Dedicated to promoting bird conservation in Uganda, ACUS collaborates with organizations like BLI to advocate for the preservation of birds, habitats, and ecological assets. Established in 2022, ACUS is the largest dedicated regional bird society in East Africa. Following their visit to BLI, ACUS recommended capacity building for site guides at the institute through the provision of birding education materials, such as guidebooks and binoculars. Given than BLI USA is headquartered in Indiana, we found a partnership with Indiana Audubon and were the grateful recipients of a 2022 Mumford & Keller Grants and Scholarship Program award. This grant facilitated the purchase of 12 sets of Gosky 10x42 Roof Prism HD binoculars for Lazarus' Forest Conservation. BLI Uganda, inspired by Pope Francis' encyclical Laudato S, addresses environmental degradation through a new integrated educational program for a lifestyle change through education and experience. Lazarus' Forest serves as a unique educational center, allowing experiential learning about conservation, integrated ecosystems, and wildlife. BLI welcomes visitors for guided tours and overnight stays, providing a chance to explore the forest's flora, fauna, and birdlife. Early this year, I personally delivered the binoculars to BLI in Uganda. What I found on this trip both surprised and delighted me. My first stop was at the home and guesthouse of one of our founders, Fr. Emmanuel Katongole. The guesthouse sits adjacent to the Entebbe Zoo, where birds are free to fly around the area and over adjoining Lake Victoria. I shared the binoculars with him and his staff to explore a bit. Most of them had never seen binoculars before, including his aunt who is well into her 80s. They were in awe of the birds, flora, and even the moon through the binoculars! It was a heart-warming experience. On to the BLI campus the next day, I delivered the binoculars to our resident Forester and Taxonomist, Chris Mukwaya, who promptly catalogued them into the Institutes inventory. The next day Chris led a class on how to use the binoculars. Our students come from the poor rural areas nearby, and struggle to pay even basic school fees. Equipment such as binoculars are non-existent. Chris led off with basic instructions on how to assemble, care and use them. He emphasized that they were important equipment in conservation, and with his wry sense of humor told the students should they break them, their father would have to sell a cow to buy new ones. The point was made the students take great care with the equipment. They were taught how to move through the forest quietly so as to not disturb the birds. And then they were off, with the assignment to find 3 common bird species, and bring back MUMFORD & KELLER HIGHLIGHT photos of birds he had not listed, using the cell phone adaptors. (Yes, the students have basic cell phones with photo capabilities, but no wifi or smartphones.) One hour later, they returned, assignments successfully completed, followed by a discussion on the varying bird species they discovered. The binoculars are integrated into the curriculum as part of the Forestry Management & Ecological Protection Cycles. In 2022, ACUS identified 100 bird species in BLI's forest, a number that has since risen to 161. The students are enjoying birdwatching and have created their own club and their theme is "their life is our life." We are in the process of building the Lazarus Forest Information Center an ecological and information center at the entrance to our 300 acre forest. From here, students and tourists check out the binoculars for a small fee to explore the forest. To date, we have had a few tourists, and word is spreading, as I even had an inquiry from a US Audubon Society member who had heard about us and wanted to visit! Interns and students from Wheaton College, the University of Notre Dame, and Duke University visited this summer and had in- depth classes in ecological and forest preservation using the binoculars. To sustain life on a healthy planet, interconnected ecosystems must be protected. Uganda, known as the "Pearl of Africa," is rich in biodiversity, ranking second richest for mammals in Africa, second for birds in Africa, and seventh for higher plant species in the world. BLI, Indiana Audubon, and similar partnerships worldwide exemplify how conservation efforts can transcend borders. We are proud to work together in appreciation of the beauty that surrounds us. We invite you to follow our progress on our website (bethanylandinstitute.org) and social media channels, as we continue restoring the ecosystem and educating our students and visitors on the care of our common home. We also invite you to visit us in Uganda, for what we promise to be an enriching experience, for both the visitor and our campus staff and students. Editor's Note: This report has been edited for length and clarity. FEB. MARCH 2024 4 FIE LD NOTES FRO M A NATU R A LIST By Joni James Nature journaling involves using words, pictures, and numbers to record your observations, connections, and experiences with nature on paper. It is not about being an artist. Observation, curiosity, and creativity are skills that you can develop. There is no wrong way to do it just have fun! Questions about this series or nature journaling? Email Joni at heronwatch0@gmail.com. 5 CARDINAL | INDIANA AUDUBON SOCIETY Collaborative Efforts R E S U LT I N AN Indiana First BY KRIS TIN S TR AT TON On June 17, 2023, The Dunes-Calumet Audubon hosted a field trip led by Julie Bonnema, with experienced birders Angie Archer and me, among the fifteen participants. It was a beautiful day for birding, filled with old and new birding friends. The day at Kankakee Sands offered perfect weather, with clear skies, a slight breeze, and temperatures in the low 70s. Julie had a few key locations in mind, each with target birds for the group to observe. Along the way, we encountered some little surprisestwo of those being at the corner of where units F, D, G and the Beaver Lake Nature Preserve meet. We stopped there to look and listen for Blue Grosbeaks. I ventured away from the group to unit G to scan and listen while Julie stayed by units F and D to do the same. Upon scanning the restored grassland prairie habitat, managed by The Nature Conservancy, Julie came across a late Clay-colored Sparrow. She called me and Angie over to see it and confirm the identification, saying, Kristin knows these well." As I was walking over, I heard the raspy buzzzz buzzzz call and confirmed that it was indeed a Clay-colored Sparrow. I snapped one photo of the bird perched with food. Julie, Angie, and I all assisted in helping the field trip attendees see the bird. We were excited because it was mid-June and Clay-colored Sparrows are migrants in Indiana, and spring migration had been over for a few weeks. As we were observing the sparrow perched with a tiny inchworm in its bill, we heard a second buzzzz buzzzz call come from the west of unit F. Two Clay-colored Sparrows! That, in itself, is a treat. We spent a considerable amount of time observing the birds and their behavior. The first bird was perched in the same spot with food, flying down to the same location four times. Julie also witnessed the bird carrying food to the same spot multiple times. Both birds were singing and doing their chip calls. At one point, Julie saw one fly away, and the other arrived with food at the same spot. Angie kept all our eBird checklists that day and logged the Clay-colored Sparrows as breeding birds. At this time, I dont think any of us had realized that we had discovered the first breeding Clay-colored Sparrows in Indiana. We needed to get the trip moving along and also didnt want to spend too much time at that spot as to disturb the sparrows. (By the way, we never heard or saw Blue Grosbeaks that day.) On our way to our next spot, we had to abruptly stop for a singing Connecticut Warbler, also a late migrant. That's a whole different story for another time. REMARKABLE RECORDS I consulted with Jeff McCoy about our observations of Connecticut Warbler and Clay-colored Sparrow, inquiring if either bird breeds in Indiana. Jeff confirmed that both species were considered late. He also indicated that breeding Clay-colored Sparrows had never been confirmed in Indiana. Upon sharing this information with Julie and Angie, we decided we had better do some more observations. The next morning, Julie and I returned to do a mini "big sit, but strong south winds altered the weather, keeping all the birds hunkered down and quiet. Despite spending four hours watching and waiting for Clay-colored Sparrow activity, we only heard one buzzzz buzzzz song and observed a single sparrow popping up from unit D, flying over and perching in a tree along the road. A few weeks later, Brad Bumgardner let us know that he had consulted with Allisyn-Marie Gillet and Ken Brock regarding our observations. He clarified that, according to the Breeding Bird Atlas protocol, a bird carrying food is considered a confirmed breeding record. Dr. Brock provided additional context, citing four other summer records in Indiana, including a pair in June 2016 at Kankakee Sands Unit D. Other records include a territorial Clay-colored Sparrow in St. Joseph County in July 1989, Pigeon River Fish and Wildlife Area in June 2011, and again in St. Joseph County in June 2013. So, on June 17, 2023, Julie, Angie, and I collectively had a hand, eye, and ear in finding the first confirmed breeding Clay-colored Sparrow in Indiana. The Clay-colored Sparrow and Connecticut Warbler sightings came at a time when I lacked confidence, especially in birding. Some great memories happened on that field trip that I will always hold near and dear to my heart. The birds, fellow trip participants, and those we consulted with served as a reminder to stay positive, be mindful, have confidence in aspects of life. This experience reignited my passion for sharing my birding knowledge with others. Caught in Action: Kristin Stratton swiftly snapped this photo of the Claycolored Sparrow carrying food, confirming the breeding sighting. FEB. MARCH 2024 6 SABBATICAL WELL-SPENT THE THRILLS OF A LOWER 48 BIG YEAR BY DAVID BENSON T heres no way we are going to find this thing. I was sure of it. The trails were packed with tourists, hiking noisily all over the mountain side, completely oblivious to the possibility of a Sooty Grouseand that there was such a thing as a Sooty Grouse. Lets just keep scanning, my friend, Maikel Wise, who joined me on this western trip to Rainier, Washington, and the Pacific coast, suggested. So, we continued scanning the heath subshrubs. Even at a distance, Ive got a decent shot at picking one out. Ive been doing research on ptarmigan in Glacier National Park, Montana, for decades and have a pretty good search image for grouse. However, my internal clock was telling me that we needed to cut our losses and move on. We spent way too much time in Seattle traffic yesterday and needed to get to the coast where the majority of my target species were. Sooty was the only species I needed in Rainier. On the other hand, it would be a lifer I cant give up too easily! Keep searching Got it! Maikel exclaimed. Wait, is that a grouse or a ptarmigan? Maikel had found the bird in the deep 7 CARDINAL | INDIANA AUDUBON SOCIETY willows, but not too far from the trail. Getting it in my binoculars, I said, Thats no ptarmy! Its a Sooty! Yes! #468 for the year, and one that I was definitely not guaranteed to get. Im a biology professor at Marian University in Indianapolis. I decided to do this big year in early March 2023 when I was granted a sabbatical for the fall. This gave me freedom for a few more fall trips than I would be able to do with my typical class schedule. My travels fit the times I had available. Spring break I went to northern Florida and south Georgia, adding 80 species to the year. Then, in early May I went to Utah to watch my middle daughter graduate from Westminster College. I went a day early so I could bird, adding 32 species. From there, I flew straight to Arizona, primarily birding in the Madera Canyon and the Portal areas for a week. There, I added 90 species and had so much fun. Everyone there is a birder, so if you are missing a target, just ask. I asked about Spotted Owl and was told the exact tree it would be sitting in during the day, along with Red-faced and Olive Warblers and Mexican Chickadees singing nearby. Returning to Indiana during the peak of spring migration and doing my required guide to go for the heron on private property south of the border wall, and he graciously invited me, along with Katey Powell and Kyle Wiktor. This was the highlight of the trip with fantastic views of some of the rarest of the rare, like the heron (#600 for the year), Brown Jays (so big that at first glance, I thought they were hawks), Morelets Seedeaters, Red-billed Pigeons, and Rose-throated Becards. Texas yielded 36 new birds for the year leaving me at 611 for the lower 48. Among the 611 species I saw in 2023, 75 were entirely new to me. Which was my favorite? For me, often its the chase more than the bird itself that endears me to it. For example, the Black Rosy-Finch is a beautiful bird on its own. But, the fact that I had to search along the Beartooth Highway, a 11,000-foot high stretch of road in the spectacular mountains just outside of Yellowstone National Park, made this species my favorite of the year. I searched and searched with mind-blowing views in all directions, and when dusk was falling and Id almost given up hope, I heard one calling outside the car window. A brief hunt located three individuals among the alpine rocks, species #446 for the year. Over the course of the year, I learned an enormous amount, working hard on tricky identifications, better understanding birds' movements, and simply grasping at a more visceral level bird ecology, which I will be able to pass on to my Marian University students. This crazy goal was the driving, fun force in my life last year, and I highly recommend it! David Benson, Ph.D., is Marian University's senior biology professor. He finished 2023 on a high note, ranking 12th for identifying an outstanding 611 bird species across America's Lower 48 states during the 2023 calendar year. Great Gray Owl (left) and White-tailed Ptarmigan (right). David Benson annual (carbon neutral) birdathon for the Amos Butler Audubon Society with Wesley Homoya resulted in an additional 50 new species. During my Montana summer I was able to bird in Glacier, of course, but also take day trips a few times around Montana, seeing lots of amazing views, habitats, and birds, and adding 70 species. Because I was already close by, I decided to go straight from Montana to Washington state in mid-September, then take a week in California, mostly in the San Francisco and Los Angeles areas, where I took a second pelagic trip. The first was off Cape Hatteras in May, adding 8 species. This one was out of Half Moon Bay and was super productive. Even though I had spent a lot of effort trying to learn the species that were most common, brown and white birds zipping by the bouncing boat at great speed take practice to identify in the field. Thankfully, pelagic tours include expert spotters on board. All told, the Washington and California trip added 73 species. In October, my wife Kara and I traveled to south Florida and the Keys, including Dry Tortugas National Park, to celebrate our 30th anniversary. Despite not being a birder, she enjoyed the thrill of adding to the burgeoning year list, bumping it from 542 to 560. Who wouldnt enjoy joining the experts at the Florida Keys Hawkwatch, seeing hundreds of peregrine falcons along with my targets, Short-tailed Hawk and Mississippi Kite? Thanksgiving found my family in North Carolina with friends. After that, I dropped off the rest of the family at the airport and drove overnight to Cape Cod, spending hours on cold, windy beaches looking, unsuccessfully, for Dovekie and Manx Shearwaters. I loved it, and I found Great Shearwater, Razorbills, and eiders, plus 10 other new species along the Massachusetts and southern Maine coasts. The last best trip of the year was to Texas, where, thanks to Jeremiah Oden, the Bare-throated TigerHeron awaited. Jeremiah had already booked the FEB.MARCH 2024 8 Q T ELL US A LIT TLE ABOUT YOURSELF. FAMILY? WHERE DO YOU LIVE AND WORK? Photo courtesy of Leah Baker I recently moved from Franklin, Indiana, to a cabin in the woods in Nashville, Indiana. I live here with my husband of almost 12 years and a senior brown tabby cat who showed up on our porch shortly after we moved. I have worked for the Indiana State Bar Association since 2017. We shifted to about 75% remote work after the pandemic, allowing me to spend most days working near my feeders and the mature trees that surround our property. I started painting birds and other fauna in 2020. Currently, I am selling some of my work and have recently been accepted to showcase it at B3 Gallery in Nashville. I hope to someday teach a class about nature journaling to inspire people to explore the natural world through creativity. Additionally, I enjoy traveling, backpack camping, and hiking. Q W HAT GOT YOU INTERES TED IN BIRDS AND BIRDING? meet a MEMBER LEAH BAKER In my early twenties, a close friend passed away unexpectedly and tragically. I remained close to his family, and his father was using some Native American rituals and lore to navigate the grieving process. Through this, he mentioned that his sons guardian was a red-tailed hawk. Naturally, I started seeing Red-tailed Hawks everywhere. This led me to actively look for them, observing their habits, and, of course, making observations of other species. Before long, my husband purchased my first birding book and a Vortex monocular. I was captivated and delighted by every new species. Eventually, after watching me take photos with my phone and monocular, my husband also gifted me a Sony a6000, and I slowly began documenting in a more intentional way. Q W HAT INSPIRED YOU TO INCORPOR ATE BIRDS INTO YOUR AR T WORK? Like many, I felt like I was losing my mind during the early days of the pandemic. Everything was shrouded in uncertainty, and it was very isolating. I needed something to occupy my time. Im not sure what possessed me to take up painting, but in short order, I enrolled in a fabulous online course with the amazing artist Liz Clayton Fuller through Cornell Lab. The topic was nature journaling. Over the next few months, I took as many courses on painting or natural journaling as I could find, and eventually, the training wheels came off. I started painting birds and nature journaling on my own. Essentially, I said, how else can I bird? And the answer was through artistic study. To me, it's another way to explore anatomy, appearance, behaviors, etc. Painting was my sourdough. Q W HAT'S YOUR FAVORITE BIRDING DES TINATION, BOTH IN INDIANA AND OUTSIDE THE S TATE? In Indiana, I like birding hotspots that are off the beaten path. I prefer land trusts because they are usually less crowded. Beanblossom Bottoms is my current favorite. I love the variety that wetlands bring. Outside of Indiana, one of my all-time favorite birding destinations is Point Lobos in California. Bird Island houses a huge colony of cormorants, and Blackcrowned Night Herons nest in the cliffs. Peregrines and Osprey soar. Oystercatchers, guillemot, Brown Pelicans, quail... it is truly extraordinary where the land meets the sea. Q D O YOU HAVE A FAVORITE BIRD OR GROUP OF BIRDS? Raptors have always been my favorite. There is something that is just mesmerizing about their piercing eyes, their athleticism, and their resilience. My current favorite bird is the Eastern Screech-Owl. I was lucky to have a nesting family visit my yard each year when I lived in Franklin. It brought me so much joy to observe them each year. 9 CARDINAL | INDIANA AUDUBON SOCIETY Resplendent Quetzal by Francesco Veronesi/Wikimedia Commons February 1-10, 2025 Join Indiana Audubon for an exciting birding adventure in Guatemala, the northernmost country in Central America. Discover the countrys rich avian diversity with nearly 750 bird species, including 35 Mesoamerican endemics. From lowland rainforests to high-elevation cloud forests, explore Guatemalas stunning landscapes featuring volcanoes, lakes, Mayan culture, and ancient ruins. This tour takes you through misty cloud forests, colonial Antigua, and offers glimpses of specialty bird species like the Horned Guan, Resplendent Quetzal, and the captivating Pink-headed Warbler. Dont miss this spectacular trip to one of Central Americas most diverse nations! El Quetzal boat by Wesley Homoya INCLUDED IN THE TOUR EXPERIENCE: IAS leader and guide coordinating trip logistics before and during the tour, LEARN MORE Lake Atitlan by Luk Janika/Pixabay Guatemala: Go Birding In the Shadow of Volcanoes including all eBird checklist submissions. Specialized Natural Selections Driver and local guides. Transportation with fuel included. All lodging for 10 days, 9 nights in listed hotels. All meals, beginning with dinner on day 1 and ending with breakfast on day 10. All entrance fees to reserves and national parks. This tour is organized by Indiana Audubon in partnership with Natural Selections Tours. Resort by Wesley Homoya $4,538.00 per person $250 Single supplement Wine-throated Hummingbird by Wesley Homoya INDIANA AUDUBON SOCIETY ADVERTISE IN THE CARDINAL For more information about The Cardinal and advertising rates, visit indianaaudubon.org/the-cardinal.  ...

... 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.  ...

... ONLINE FIRST This is a provisional PDF only. Copyedited and fully formatted version will be made available soon. ISSN: 0015-5659 e-ISSN: 1644-3284 A unique variation of the jugular veins and its clinical significance Authors: McKenzie Young, Alexis Zavitsky, Jillian Niceley, Gabriella Battiston, Sumathilatha Sakthi-Velavan DOI: 10.5603/fm.98618 Article type: Case report Submitted: 2023-12-23 Accepted: 2024-02-13 Published online: 2024-02-23 This article has been peer reviewed and published immediately upon acceptance. It is an open access article, which means that it can be downloaded, printed, and distributed freely, provided the work is properly cited. Articles in "Folia Morphologica" are listed in PubMed. Powered by TCPDF (www.tcpdf.org) CASE REPORT A unique variation of the jugular veins and its clinical significance McKenzie Young et al., Variant veins of the neck McKenzie Young, Alexis Zavitsky, Jillian Niceley, Gabriella Battiston, Sumathilatha SakthiVelavan Division of Biomedical Sciences, Marian University College of Osteopathic Medicine, Indianapolis, IN, 46222, USA Address for orrespondence: Sumathilatha Sakthi-Velavan, Associate Professor of Anatomy, Division of Biomedical Sciences, Marian University College of Osteopathic Medicine, 3200 Cold Spring Rd, Indianapolis, IN 46222, USA; e-mail: ssakthivelavan@marian.edu, tel./fax : +1(317) 955-6245 ABSTRACT The authors report a rare variation of the anterior jugular and internal jugular veins in a 78-year-old male donor. An enlarged and curved left anterior jugular vein (AJV) was formed as the continuation of the left common facial vein (CFV). The left AJVs diameter was wider than the internal jugular vein (IJV) and measured around 5 mm greater than the IJVs diameter and a channel connected the two veins. The right AJV and CFV continued from the two divisions of the right facial vein. The right AJVs diameter was smaller than the right IJVs diameter. The right external jugular vein was absent. No concurrent pathology supported the abnormal dimension of the left AJV and the findings were indicative of a variant anatomy. These variations have rarely been reported and have important clinical correlations. Failed IJV cannulation may result if the variant neck veins are missed. However, variant veins may serve as collateral channels and patch material in IJV reconstruction, carotid angioplasty, and ventricular-jugular shunts. Keywords: jugular veins, anatomical variation, catheterization, venous anastomoses, communicating vein, bilateral variation INTRODUCTION 1 The veins of the head and neck drain the brain, face, and neck [1]. The anterior jugular vein (AJV), formed by the confluence of submandibular veins near the midline, drains into the external jugular vein (EJV) or subclavian vein (SCV) [2]. The internal jugular vein (IJV), the continuation of sigmoid sinus, terminates by uniting with the SCV to form the brachiocephalic vein [2]. The facial vein (FV) unites with the anterior division of the retromandibular vein (RMV) to form the common facial vein (CFV) which drains into the IJV [2]. The posterior division of the RMV unites with the posterior auricular vein to form the EJV which drains into the SCV [2]. Anatomical variations of neck veins are common because of the regions robust vascular network [3]. A clear understanding of the venous variations helps prevent complications like pneumothorax, muscular damage, and carotid artery puncture during procedures [1, 4]. Commonly reported variations in the superficial neck veins include the EJV (most variable), posterior EJV, AJV, and the thyroid veins [1]. Studies reported 26% of 104 patients had either a unilateral or bilateral IJV variation, 0.3 to 2.0% of 330 patients had facial vein (FV) variations, and 6 of 35 (17.1%) donors had CFV variations [3, 5, 6]. This case report explores a unique course and drainage of the AJV and CFV bilaterally, and their clinical significance. CASE REPORT An atypical variation in the formation, dimensions, and the drainage of the left anterior jugular vein was observed in a 78-year-old male donor during routine dissection in the anatomy laboratory. Detailed dissection preserving the veins and related neurovasculature was implemented. The origin, drainage, and tributaries of other neck veins, were explored. The IJV was traced to the cranial cavity and the intracranial sinuses were examined. The enlarged neck veins were opened and the lumina were examined for thrombosis or other pathology. The veins were painted blue with acrylic paint (Apple Barrel Matte Bright Blue Acrylic Craft Paint) and photographed. The vein circumferences were measured using digital calipers at the level of the thyroid cartilage, and diameters were calculated using the equation, circumference/ (Table 1). The left AJV was found to be greatly enlarged and curved on the anterolateral aspect of the neck. It was formed as the continuation of the left CFV and a communicating vein connected it to the left IJV (Figs. 1, 3B). The left AJV united with the left EJV before draining into the SCV at the level of the clavicle. The left IJV was smaller than the left AJV and its counterpart on the right. A communicating vein connected the left AJV and IJV. The left EJV was the continuation of the posterior division of the RMV, while the anterior division of RMV joined the facial vein to form the CFV. The right jugular system differed from the left. The right AJV appeared normal and measured smaller than the left AJV. The right FV terminated as two divisions (Figs. 2, 3A). The posterior division continued as the right CFV, while the anterior division continued as the AJV. Like the left 2 side, a communicating channel connected the AJV and IJV. The right IJV appeared dilated. The right EJV was absent. The lumina of neck veins showed no pathology. DISCUSSION The head and neck veins begin developing in the embryo during week 4 from the cardinal venous system [2]. The ventral pharyngeal vein (VPV), which drains the mandibular and hyoid arches, is the first recognizable vein [2]. The VPV drains into the anterior cardinal vein, which becomes the IJV [2]. Several venous channels from face and neck drains into the IJV [2]. A capillary plexus from the neck forms the EJV [2]. The convergence of the superficial submandibular veins forms the AJVs, which empty into the EJVs [2]. The variant anatomy noted in this case may be due to the lack of or malformation of the venous anastomotic channels. Nayak reported a case in which the FV continued as the AJV in the presence of a communicating vein between the AJV and IJV, unilaterally [7]. Kumar and Baidya reported a case in which the CFV continued as the AJV with no communicating vein between the IJV and AJV, unilaterally [8]. However, this case is unique since bilateral variations with the left CFV and the anterior division of the right FV continuing as the left and right AJV respectively were present. The dimensions of the IJVs and left AJV in this case were noteworthy. Tartire et al. found that the right IJV was larger than the left in 75 to 80% of 190 adult patients [9]. The mean diameter of the right IJV was 17 5 mm, and the mean diameter of the left IJV was 14 5 mm at the level of the cricoid cartilage [9]. In this case, the right IJV was larger than the left IJV, supporting the study. However, the right IJV had a larger than normal diameter above the cricoid cartilage, and the left IJV had a smaller than normal diameter (Table 1). Hojaij et al. reported that the AJVs diameter was 5 mm in 40% of 30 cadavers studied; however, the AJV is typically smaller than the other jugular veins [2, 10]. In this case, the left AJV was larger than the other jugular veins on the left (Table 1) which was likely due to its communication with the other veins. Enlarged superficial veins may indicate a pathological process [7, 11]. Additionally, if the AJV is larger than the EJV, then thrombosis of the EJV may be present [7, 11]. However, no evidence of thrombosis was found in this case, which is suggestive of a variant anatomy [11]. It is likely that the enlarged right IJV and left AJV formed secondary to the smaller left IJV to provide adequate drainage of the head and neck. An absent EJV, as noted on the right side of this donor, may be due to the lack of or malformation of anastomotic channels, such as that between the cephalic vein and FV [2, 12]. In the case of an absent EJV, the veins that typically form the EJV open into the IJV as noted in this case [12]. Catheter mispositioning or failed IJV cannulation may result if the variant superficial veins are not noted [11]. Utilizing ultrasound may help differentiate the IJV from a prominent AJV during 3 cannulation, to prevent unintentional injury [11]. A discernible AJV may serve as a collateral pathway for intracranial venous drainage [11]. Specifically, in the case of unilateral occlusion of brachiocephalic vein and when performing IJV reconstruction, an enlarged AJV is considered useful [11, 13]. Normally, the union of the CFV and IJV is just superior to the carotid bifurcation, providing the surgeon with a consistent landmark [4]. Adequate and accurate access to the bifurcation is essential because it is the most common location of atherosclerosis [14]. When the left CFV does not unite with the IJV, as noted in this case, it cannot be used as an accurate landmark. During carotid angioplasty and endarterectomy, the anterior border of IJV may be helpful in locating the incision site to open the carotid sheath [4]. While an enlarged IJV helps in an easier localization of incision site, it poses an increased risk of injury to the IJV compared to a procedure with a normal or smaller IJV. During cannulation, the right IJV is preferred because it provides direct access to the superior vena cava, is more commonly the dominant hand side, has a lower complication rate, and is typically larger than the left IJV [9]. The larger than normal right IJV, as noted in this donor, provides easier access for cannulation. CONCLUSIONS The presence of variant head and neck veins, reported in this case, can complicate patient care while also providing alternatives if other pathologies or variations are not noted. Without knowledge of a patients variant anatomy, the physician may not be successful in cannulation, catheterization, carotid angioplasty and endarterectomy, and other procedures. Article information and declarations Ethics statement The Marian University Institutional Review Board declared that the study did not need review or approval and cleared the study since the research was on a cadaver (IRB#B23.109). Author contributions McKenzie Young: Study conception, design, dissection, data collection, manuscript preparation and review and preparation of figures 1, 2, and 3. Alexis Zavitsky: Study conception, design, dissection, data collection, manuscript preparation and review and preparation of table 1. Jillian Niceley: Study conception, design, manuscript preparation and review and preparation of Table 1. Gabriella Battiston: Study conception, design, manuscript preparation and review. 4 Sumathilatha Sakthi-Velavan: Conceptualization, guidance, reviewing, editing and funding acquisition. Acknowledgments The authors would like to thank the donor and his family for their generosity so that anatomical research could be performed. The authors also thank Marian University College of Osteopathic Medicine for the support in conducting the research. Conflict of interest The authors declare no conflict of interest. Funding The study was supported by the Marian University Research & Scholarship Administration. REFERENCES 1. Dalip D, Iwanaga J, Loukas M, et al. Review of the variations of the superficial veins of the neck. Cureus. 2018; 10(6): e2826, doi: 10.7759/cureus.2826, indexed in Pubmed: 30131919. 2. Standring S. Development of the head and neck. In: Standring S. ed. Gray's anatomy: the anatomical basis of clinical practice. Elsevier Churchill Livingstone, Edinburgh 2021: 273 291. 3. Bertha A, Rabi S. Anatomical variations in termination of common facial vein. J Clin Diagn Res. 2011; 5(1): 2427. 4. Gupta V, Tuli A, Choudhry R, et al. Facial vein draining into external jugular vein in humans: its variations, phylogenetic retention and clinical relevance. Surg Radiol Anat. 2003; 25(1): 3641, doi: 10.1007/s00276-002-0080-z, indexed in Pubmed: 12819948. 5. Bondaz M, Ricard AS, Majoufre-Lefebvre C, et al. Facial vein variation: implication for facial transplantation. Plast Reconstr Surg Glob Open. 2014; 2(7): e183, doi: 10.1097/GOX.0000000000000134, indexed in Pubmed: 25426366. 6. Lin BS, Kong CW, Tarng DC, et al. Anatomical variation of the internal jugular vein and its impact on temporary haemodialysis vascular access: an ultrasonographic survey in uraemic patients. Nephrol Dial Transplant. 1998; 13(1): 134138, doi: 10.1093/ndt/13.1.134, indexed in Pubmed: 9481729. 7. Nayak BS. Surgically important variations of the jugular veins. Clin Anat. 2006; 19(6): 544 546, doi: 10.1002/ca.20268, indexed in Pubmed: 16372344. 8. Kumar S, Baidya R. Termination of the common facial vein into the anterior jugular vein: a rare case report. Int J Anat Res. 2018; 6(3.2): 55015503, doi: 10.16965/ijar.2018.255. 9. Tartire D, Seguin P, Juhel C, et al. Estimation of the diameter and cross-sectional area of the internal jugular veins in adult patients. Crit Care. 2009; 13(6): R197, doi: 10.1186/cc8200, indexed in Pubmed: 20003190. 5 10. Hojaij FC, Santos L, Moyses R, et al. Anatomy of the anterior jugular veins: anatomical study of 30 cadavers. Arch Head Neck Surg. 2022; 51: e20220005, doi: 10.4322/ahns.2022.0005. 11. Schummer W, Schummer C, Bredle D, et al. The anterior jugular venous system: variability and clinical impact. Anesth Analg. 2004; 99(6): 16251629, doi: 10.1213/01.ANE.0000138038.33738.32, indexed in Pubmed: 15562044. 12. Cvetko E. A case of left-sided absence and right-sided fenestration of the external jugular vein and a review of the literature. Surg Radiol Anat. 2015; 37(7): 883886, doi: 10.1007/s00276-014-1398-z, indexed in Pubmed: 25432662. 13. Katsuno S, Ishiyama T, Nezu K, et al. Three types of internal jugular vein reconstruction in bilateral radical neck dissection. Laryngoscope. 2000; 110(9): 15781580, doi: 10.1097/00005537-200009000-00034, indexed in Pubmed: 10983966. 14. De Syo D, Franji BD, Lovricevi I, et al. Carotid bifurcation position and branching angle in patients with atherosclerotic carotid disease. Coll Antropol. 2005; 29(2): 627632, indexed in Pubmed: 16417173. 6 Figure 1. Left head and neck venous system, 1 thyroid cartilage, 2 internal jugular vein, 3 external jugular vein, 4 posterior division of retromandibular vein, 5 retromandibular vein, 6 anterior division of retromandibular vein, 7 facial vein, 8 common facial vein, 9 anterior jugular vein, 10 subclavian vein, 11 communicating vein between IJV and AJV, 12 common trunk of superior and middle thyroid veins, 13 brachiocephalic vein, 14 common carotid artery, 15 submandibular gland. 7 Figure 2. Right head and neck venous system, 1 thyroid cartilage, 2 internal jugular vein, 3 retromandibular vein, 4 facial vein, 5 common facial vein, 6 anterior jugular vein, 7 communicating vein between IJV and AJV, 8 common carotid artery, 9 submandibular gland. 8 Figure 3. A. Schematic of right head and neck venous system, B. Schematic of left head and neck venous system. 1 submandibular gland, 2 omohyoid, 3 sternocleidomastoid, 4 internal jugular vein, 5 external jugular vein, 6 posterior division of retromandibular vein, 7 retromandibular vein, 8 anterior division of retromandibular vein, 9 facial vein, 10 common facial vein, 11 anterior jugular vein, 12 subclavian vein, 13 communicating vein between IJV and AJV, 14 common trunk of superior and middle thyroid veins, 15 brachiocephalic vein, 16 common carotid artery, 17 posterior division of facial vein, 18 anterior division of facial vein. 9 Table 1. Jugular vein measurements. Diameter at the level of Diameter at the level of angle Left AJV IJV EJV Righta AJV IJV thyroid cartilage of mandible 8.91 mm 3.82 mm 5.10 mm 2.55 mm 1.91 mm 15.92 mm 11.46 mm a) No EJV was observed on the right side. AJV anterior jugular vein; EJV external jugular vein; IVV internal jugular vein 10  ...

... 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.  ...

... Editorial 17 January 2024 DOI 10.3389/fendo.2024.1347765 TYPE PUBLISHED OPEN ACCESS EDITED AND REVIEWED BY Ralf Jockers, Universite Paris Cite, France *CORRESPONDENCE Michela Rossi michela1.rossi@opbg.net 01 December 2023 08 January 2024 PUBLISHED 17 January 2024 RECEIVED ACCEPTED CITATION Rossi M, Lowery JW and Del Fattore A (2024) Editorial: Genetic and molecular determinants in bone health and diseases. Front. Endocrinol. 15:1347765. doi: 10.3389/fendo.2024.1347765 Editorial: Genetic and molecular determinants in bone health and diseases Michela Rossi 1*, Jonathan W. Lowery 2,3,4,5,6 and Andrea Del Fattore 1 1 Bone Physiopathology Research Unit, Translational Pediatrics and Clinical Genetics Research Division, Bambino Ges Childrens Hospital, IRCCS, Rome, Italy, 2 Division of Academic Affairs, Marian University, Indianapolis, IN, United States, 3 Department of Physiology & Pharmacology, College of Osteopathic Medicine, Marian University, Indianapolis, IN, United States, 4 Bone & Muscle Research Group, Marian University, Indianapolis, IN, United States, 5 Indiana Biosciences Research Institute, Indianapolis, IN, United States, 6 Indiana Center for Musculoskeletal Health, Indiana University School of Medicine, Indianapolis, IN, United States KEYWORDS bone, osteoclast, osteoblast, bone disease, gene COPYRIGHT 2024 Rossi, Lowery and Del Fattore. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms. Frontiers in Endocrinology Editorial on the Research Topic Genetic and molecular determinants in bone health and diseases Alterations of bone remodeling impact skeletal integrity lead to excessive or impaired bone resorption as well as reduced or disorganized bone formation (1, 2). This Research Topic focuses on the identication of genetic and molecular determinants involved in both bone health and diseases. In this editorial, we highlight studies on rare diseases presented in the Research Topic with the aim of better understanding their etiopathogenesis and opening the way for the identication of new therapeutic approaches. Xiang and Zhong summarized the recent studies regarding the molecular and cellular mechanisms leading to the progressive osteolysis and angiomatous proliferation in Gorham-Stout disease (GSD), which is a very rare disease that is also known as Vanishing Bone Disease. GSD is characterized by severe osteolytic bone destruction but lacks specic diagnostic markers and therapy (3). The information presented by Xiang and Zhong provides an important update on the condition and presents ideas for new therapeutic approaches for this rare disease. Cinque et al. published an elegant study on hypophosphatasia (HPP), a rare genetic disease affecting bone and teeth mineralization with multisystemic manifestations involving the nervous system, musculoskeletal apparatus, and kidneys, due to ALPL mutations. The authors reported the genetic analysis performed on 33 patients, identifying eight novel variants of ALPL gene. These results associated with the detailed clinical description increase the knowledge of this rare condition. Osteogenesis imperfecta (OI), also known as brittle bone disease, is a clinically and genetically heterogeneous disorder of connective tissue and is identied by bone dysplasia and fragility (4). In this Research Topic, Paduano et al. reported the results obtained by next-generation sequencing (NGS) analysis of 10 patients, comprising 7 male and 3 female patients from 7 families, all from the Puglia Region in South Italy. The authors identied novel rare pathogenic variants in type I collagen-encoding genes (COL1A1 and COL1A2). 01 frontiersin.org Rossi et al. 10.3389/fendo.2024.1347765 10 years of age then plateaus until old age, with a trend of bone turnover markers similar to that of humans. In conclusion, the papers published in this Research Topic underline how investigating bone diseases and animal models represent a way to nd new determinants of bone physiology and also allow the identication of new therapeutic approaches. In another study regarding OI, Lim et al. described the effects of a missense variant of MBTPS2 which encodes the site-2 protease, a Golgi transmembrane protein that activates membrane-tethered transcription factors in aborted male fetus with micromelia particularly of the lower limbs, a narrow thorax, and defective ossication of calvarium. The authors performed in vitro studies on mutated MBTPS2 primary broblasts and found perturbations in fatty acid metabolism and collagen production. Sundqvist et al. report a case study on rare, chronic nonbacterial osteomyelitis (CNO). They described a female patient with CNO with systemic inammation, advanced malnutrition and complete deciency of myeloperoxidase (MPO). The authors reported that, although the patient did not nd benecial effects after treatment with nonsteroidal anti-inammatory drugs, corticosteroids, bisphosphonates or IL1-receptor antagonists (anakinra), the administration of TNFa blockade (adalimumab) resulted in instant resolution of the inammatory symptoms suggesting that the disease was TNFa-driven. Bone tissue is tightly connected with other organs to regulate whole physiology (5). In this Research Topic the interplay bone-liver has been reported. Huang et al. investigated whether serum liver enzymes are causally associated with bone and joint-related diseases using Mendelian randomization (MR) designs. Indeed, the positive causality between ALP and the risk of osteoporosis and rheumatoid arthritis was indicated. Moreover, the authors reported that higher levels of alanine transaminase (ALT) were associated with the risk of hip and knee osteoarthritis while no causal relationship between GGT and bone and joint-related diseases was revealed. Moreover, two further papers reported new advances in bone remodelling, using animal models. Verlinden et al. investigated how neuropilin 2 (NRP2) in osteoblasts regulates trabecular bone mass in male mice. NRP2 is a non-tyrosine kinase transmembrane glycoprotein receptor. The authors generated two different genetic models lacking Nrp2 expression in osteoblasts or osteoclasts to identify its role in the bone remodelling activity. Although loss of Nrp2 in the osteoclast lineage did not result in a bone phenotype, loss of Nrp2 in osteoblast precursors and mature osteoblasts leads to reduced cortical crosssectional tissue area and lower trabecular bone content in male mice. Li et al. performed the evaluation of bone turnover markers and DEXA (Dual-Energy X-Ray Absorptiometry) analysis in cynomolgus monkeys at different ages to establish an animal model for age-related osteoporosis in non-human primates. The authors nd that, in cynomolgus monkeys, peak BMD occurs at age Author contributions MR: Writing original draft, Writing review & editing. JL: Writing original draft. AD: Writing original draft, Writing review & editing. Funding The author(s) declare nancial support was received for the research, authorship, and/or publication of this article. MR is supported by the Fondazione Umberto Veronesi. This work was also supported by the Italian Ministry of Health with the Current Research funds. Conict of interest The authors declare that the research was conducted in the absence of any commercial or nancial relationships that could be construed as a potential conict of interest. The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the nal decision. Publishers note All claims expressed in this article are solely those of the authors and do not necessarily represent those of their afliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher. References 1. Unnanuntana A, Rebolledo BJ, Khair MM, DiCarlo EF, Lane JM. Diseases affecting bone quality: beyond osteoporosis. Clin Orthop Relat Res (2011) 469 (8):2194206. doi: 10.1007/s11999-010-1694-9 3. Rossi M, Buonuomo PS, Battafarano G, Conforti A, Mariani E, Algeri M, et al. Dissecting the mechanisms of bone loss in Gorham-Stout disease. Bone. (2020) 130:115068. doi: 10.1016/j.bone.2019.115068 2. Feng X, McDonald JM. Disorders of bone remodeling. Annu Rev Pathol (2011) 6:12145. doi: 10.1146/annurev-pathol-011110-130203 4. Martin E, Shapiro JR. Osteogenesis imperfecta:epidemiology and pathophysiology. Curr Osteoporos Rep (2007) 5(3):917. doi: 10.1007/s11914-007-0023-z 5. Yuan W, Song C. Crosstalk between bone and other organs. Med Rev (Berl). (2022) 2(4):33148. doi: 10.1515/mr-2022-0018 Frontiers in Endocrinology 02 frontiersin.org  ...