iScience ll OPEN ACCESS Article Insular giant leporid matured later than predicted by scaling Meike Köhler, Carmen Nacarino- Meneses, Josep Quintana Cardona, Walter Arnold, Gabrielle Stalder, Franz Suchentrunk, Salvador Moyà- Solà meike.kohler@icp.cat Highlights Island giants grow longer, not faster, than their continental ancestors Abundance of resources is not a prerequisite for insular gigantism In vivo labeling is fundamental for correlating histology with developmental stages Life history theory provides a framework for uncovering evolutionary mechanisms Köhler et al., iScience 26, 107654 September 15, 2023 ª 2023 The Author(s). https://doi.org/10.1016/ j.isci.2023.107654 iScience ll OPEN ACCESS Article Insular giant leporid matured later than predicted by scaling Meike Köhler,1,2,3,5,* Carmen Nacarino-Meneses,2 Josep Quintana Cardona,2 Walter Arnold,4 Gabrielle Stalder,4 Franz Suchentrunk,4 and Salvador Moyà-Solà1,2,3 SUMMARY The island syndrome describes morphological, behavioral, and life history traits that evolve in parallel in endemic insular organisms. A basic axiom of the island syndrome is that insular endemics slow down their pace of life. Although this is already confirmed for insular dwarfs, a slow life history in giants may not be adaptive, but merely a consequence of increasing body size. We tested this question in the fossil insular giant leporid Nuralagus rex. Using bone histology, we con- structed both a continental extant taxon model derived from experimentally fluorochrome-labeled Lepus europaeus to calibrate life history events, and a growth model for the insular taxon. N. rex grew extremely slowly and delayed maturity well beyond predictions from continental phylogenetically corrected scaling models. Our results support the life history axiom of the island syndrome as generality for insular mam- mals, regardless of whether they have evolved into dwarfs or giants. INTRODUCTION In their classic comparative study of island-mainland populations of giant rodents, Adler & Levins1 described under the term ‘‘island syn- drome’’ a set of covarying demographic, reproductive, behavioral, and morphological traits that they considered to be typically associated with the ecological conditions of insularity. Among these, evolution of dwarfs and giants (also called ‘‘Island Rule’’2,3), reduction in reproduc- tive output, delay in the age at sexual maturity (ASM), and a long lifespan are the most conspicuous adaptations, with a shift toward the slow end of the fast-slow life history continuum as the central axiom. The island syndrome model is supported by life history theory4 and life history models for age and size at maturity.5–7 At about the same time that life history theory was gaining ground, however, the notion that life history traits are primarily size-depen- dent8–10 and thus non-adaptive per se was spreading, leading increasingly to the use of allometric scaling approaches. Studies from this period concerned with the evolution of insular dwarfs and giants (the Island Rule2,3) therefore generally assumed that shifts in life history traits are triggered primarily by changes in body size.11–16 This led to the use of scaling models for insular dwarfs to infer shifts toward an accelerated life history with shorter growth periods,11,12,14–16 early onset of reproduction,12,13 and short lifespans.11,12,14–16 However, evidence is accumu- lating from bone and dental histology that insular dwarf taxa (elephants, bovids, cervids) from pristine (fossil) insular environments consistently delay age at maturity and extend longevity, evolving slow and not fast life histories.17–25 To explain life history evolution in insular giants, the island syndrome model and general scaling models are not mutually exclusive, since both agree that giants evolve toward the slow end of the fast-slow life history continuum. Indeed, empirical studies increasingly find adaptive trends toward life history deceleration, specifically delay in ASM and increased lifespan.26–29 However, despite decades of research, these approaches remain silent on the critical question of whether this life history shift in insular giants also evolves independently and decoupled from body mass scaling, as is the case in insular dwarfs. Among insular mammals, there is hardly a more suitable taxon than the fossil giant leporid Nuralagus rex from Menorca to resolve the question whether a shift toward a slow life history in insular giants is selected independently from body size, and thus forms part of the life history axiom of the island syndrome. Remains of N. rex come from the Late Neogene (Lower Pliocene) of Menorca,30 the 700 km2 large northernmost island of the Balearic archipelago (Spain; Figure 1). The accompanying vertebrate fauna with typical insular elements exhibits a high degree of endemism; it consists of three non-competing mammals (the giant leporid, a giant glirid, and a bat), several reptiles (one turtle, snakes, lizards, and one gecko), one amphibian, and several birds, among them two endemic owls and a corvid.31 1ICREA Pg. Lluı́s Companys 23, 08010 Barcelona, Spain 2ICP Institut Català de Paleontologia Miquel Crusafont, Edifici Z, Universitat Autònoma de Barcelona, C/ de Les Columnes, s/n., 08193 Bellaterra, Barcelona, Spain 3BABVE (Departament de Biologia Animal i d’Ecologia) Universitat Autònoma de Barcelona, 08193 Cerdanyola, Spain 4Research Institute of Wildlife Ecology, Department of Interdisciplinary Life Sciences, University of Veterinary Medicine, Savoyenstraße 1, Vienna A-1160, Austria 5Lead contact *Correspondence: meike.kohler@icp.cat https://doi.org/10.1016/j.isci.2023.107654 iScience 26, 107654, September 15, 2023 ª 2023 The Author(s). 1 This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). ll iScience OPEN ACCESS Article Figure 1. Geographic location of the Lower Pliocene fossil site Punta Nati (Menorca, Balearic Islands) The site is composed of several Late Pliocene fissure fillings, embedded in Jurassic dolomitic limestone, Miocene calcarenites, and reef deposits.30 (A) Situation of Menorca in the Mediterranean Sea (arrow). (B) Situation of Punta Nati (arrow) close to Ciutadella, Menorca; White areas: Jurassic (Dolostones); gray areas: marine Upper Miocene; The site are fissure fillings within a very narrow coastal area of the Jurassic limestone. (C) Bone breccia with remains of Nuralagus from Punta Nati (Menorca, Balearic Islands) during the preparation process. With a body mass of 8 kg,32 Nuralagus is the largest wild leporid ever known. It had a flat and small brain case, reduced orbits, and small tympanic bullae suggesting decreased visual acuity and reduced hearing capacity; a rather stiff vertebrate column, small pulmonary capacity, and short-limbed palmigrady indicate low-gear locomotion.30 These morphological and neurological adaptations are known to evolve under selection for extreme energy savings in low-predator or predator-free ecosystems where populations periodically reach high densities near carrying capacity; these traits are typically found to varying degrees in endemic island mammals such as the fossil bovid Myotragus balearicus (Mallorca, Balearic Islands;18,19,33–35), Candiacervus (Crete;36,37), and in extant though less extreme cases such as Odocoileus hemionus,21 Cer- vus nippon,22 or Amami rabbit Pentalagus furnessi.38 Pentalagus, though not a giant, is of special interest here as it is the only extant leporid with insular evolution at geological timescale. It does not only show a series of morphological characteristics similar to Nuralagus30 but it also shifted toward a slow life history. Pentalagus is reported to have the lowest fecundity rate among extant lagomorph species, with an extremely small litter size (usually 1 offspring), few breeding seasons (generally 1–2 litters/year), and late weaning.38 RESULTS By analyzing the bone histology of the insular giant leporid N. rex, we investigated whether the pace of life in insular giants is slower than predicted by continental scaling models. We used ASM and growth rate (k) as proxies of Nuralagus’ position on the fast-slow axis (or more specifically on the ‘‘reproductive-timing axis’’ sensu Bielby39). The transition from somatic growth to reproduction is a critical trade-off in mammalian life cycles.40 Several authors suggest that it is re- corded in transversal sections of long bones by a more or less abrupt change from early fast-growing tissue to extremely slow-growing tissue (external fundamental system EFS) at onset of sexual maturity.19,41–43 Hence, ASM in fossil taxa can be inferred by counting the annual lines of arrested growth (skeletochronology) deposited before formation of the EFS,41–43 by estimating the age at which the state of epiphyseal fusion in the fossil taxon is comparable to that of phylogenetic close species at onset of maturity,20 or by estimating ASM from growth functions.20 Skeletochronology is fundamental for both aging the morphological traits that record the shift from growth to reproduction (EFS and epiphyseal fusion) and for constructing the growth functions. A pre-requisite for reliable skeletochronology in fossil taxa is the correlation of bone tissue transitions with life history events in a closely related extant species with known age of sexual maturity. For this purpose, we used a large sample of Lepus europaeus that we subdivided into two groups: (i) A fluorochrome-labeled subsample (FL) and (ii) an un- labeled subsample (UnL) (Table 1; See STAR Methods for further information). 2 iScience 26, 107654, September 15, 2023 iScience ll Article OPEN ACCESS Figure 2. Bone tissue transitions and onset of sexual maturity in labeled L. europaeus (A) 1-year-old female (IPS93663a-1 ID16). Birth 31.5.2015, death 25.5.2016. Preweaning/postweaning transition from fibrolamellar complex to lamellar-primary osteon tissue at day 28 (individual without weaning LAG); Puberty: transition to lamellar tissue with simple vascularization at around 5 months of age (beginning of November); Formation of EFS with first winter LAG shortly after 8 months of age (beginning of February). left: transmitted light with labels; right: polarized light. (B) late born, 11-months-old male (IPS88859 ID15). Birth 19.7.2015, death 22.6.2016. Preweaning/postweaning transition from woven bone (FLC with some isolated Haversian Systems) to predominantly lamellar bone with primary, mostly longitudinal osteons (LPO) at day 27 (individual without weaning LAG); Puberty: transition to lamellar tissue with simple vascularization (LSV) at around 5 months of age (December); Formation of EFS with first winter LAG after 8 months of age (beginning of March). Bone tissue transitions and onset of sexual maturity Labeled L. europaeus Preweaning tissue in L. europaeus consists of a fast-growing fibrolamellar complex (FLC). Weaning is characterized by the appearance of a partial resting mark (LAG or annulus) deposited within a single day. Shortly after weaning, the tissue changes from FLC to slower growing lamellar bone with primary osteons (LPO) formed by longitudinal vascular canals (Figures 2 A, B). At puberty, bone tissue subtly shifts toward even slower growing lamellar tissue with simple vascular canals (LSV) and flattened osteocytes. Most individuals, except for early born males, reach puberty in the first half of their first winter. Sexual maturity begins after a short period of puberty. The coincidence of the onset of sexual maturity with winter depression (extremely slow growth and formation of the first winter LAG, generally between end of January and beginning of March) initiates the iScience 26, 107654, September 15, 2023 3 ll iScience OPEN ACCESS Article Figure 3. Bone tissue transitions in N. rex < 2 years old (IPS15375) (A) Partial section of the slide showing the location of the microscopic section below (red rectangle). Green arrow: weaning LAG; red arrows: 1st and 2nd(?) winter. The outer cortical area broke away (above weaning transition) deleting much of the postweaning tissue (before the first winter LAG) and the first year (between first and second(?) LAG) due to the use of aggressive acid treatment; the application of acid was necessary to free the bones from the hard cave limestone (see STAR Methods). (B) magnified micrograph of the section. Weaning transition from LPO to LSV tissue with a LAG. 1st to 2nd winter shows a decrease in the number of simple vascular channels. Area around the medullary cavity with invading Haversian systems. external fundamental system (EFS; see definition derived from our analysis of labeled European brown hares in STAR Methods). Adult in- dividuals form subsequent winter LAGs between progressively shorter annual periods of residual growth within the lamellar non-vascular (LNV) bone of the EFS. N. rex Preweaning bone tissue in Nuralagus is of LPO, seldomly (only locally) LSV type (Figure 3). All individuals show Haversian systems intruding from the medullary cavity until the weaning transition. Transition from pre- to post- weaning is recorded under polarized light as a light band of lamellar bone (annulus) that may change locally into a LAG (Figures 3 and 4D,E). The weaning discordance in all femora was deposited before the first winter LAG, clearly suggesting that Nuralagus was already weaned in the year of birth. The tissue transition from LPO/LSV to LNV tissue before the first winter growth arrest (Figure 3) shows how growth rate pro- gressively slows down from weaning onward. After the first winter LAG, there are no more primary osteons but only sparse simple vascular canals. These decrease in number until the end of the third year (Figure 4). We interpret the complete absence of vascularization (LNV tissue) as the beginning of sexual maturity (EFS; see STAR Methods). LNV tissue represents a spectrum of lowest bone growth activity, and only differences in cross-sectional tissue area between successive LAGs provide hints of the annual magnitude of bone apposition; the ESF is the last growth stage before complete arrest (no further LAGs). Compared to tissue areas as measured here (STAR Methods), radial distances between LAGs are deceptive. With growing bone diameter, bone area is distributed over an increasingly larger surface. Therefore, zones appear progressively smaller when measured linearly as radial thickness even if their area size does not change; that is, a decrease in distance between LAGs is always more accentuated than the actual decrease in bone area (STAR Methods). Although commonly used in bone histology,44 we therefore propose replacing the traditional radial 4 iScience 26, 107654, September 15, 2023 iScience ll Article OPEN ACCESS Figure 4. Bone tissue transitions and annual LAGs in L. europaeus female (IPS67842) and in N. rex female (IPS15574) (A, C, and E) Lepus female 5 years old. (B, D, and F) Nuralagus female. (E, F) entire femoral cross sections; (A, C, B, and D) magnifications. Lepus only shows two LAGs of the five expected from known age; Nuralagus shows locally a fifth LAG that is better visible under polarized light; however, it is difficult to trace at this magnification (see B). After the third LAG, there is an external fundamental system (EFS) locally traceable by 2 closely spaced LAGs in the outer cortex. Thus, the bone tissue (vascularization) and the number of annual LAGs in this individual (5), identified as a female by its size and by von Bertalanffy growth curves, suggest the beginning of an EFS at the age of roughly 3 years, which agrees with the age at maturity for females derived from their growth curves. Green arrow heads: weaning transition; red arrow heads: annual LAGs; white stars on cross sections: origin of magnifications; LAGs on cross sections: green weaning, white annual LAGs. A, B, E, F: polarized light; C, D: transmitted light. (length) measurements with area measurements since these are not subject to local variations in tissue growth rates and they are biomechan- ically more meaningful. Also in Nuralagus, the distances between the LAGs gradually decrease, making it difficult to visually determine the onset of puberty in an already slowly growing tissue. At around 3–5 years of age, locally closely spaced LAGs could indicate the onset of EFS and hence sexual matu- rity (see Figures 3 and 4; STAR Methods ‘‘External fundamental system’’). Due to the inherent optical illusion, the visual determination of ASM by LAG distances needs to be validated by growth models; these should be based on measurements of annual cross-sectional area incre- ments on transverse slides, delimited by consecutive annual LAGs (see STAR Methods). Estimating ASM from growth functions The importance of extant taxon models for inferring ASM in fossils Estimation of ASM from growth functions is a method frequently used in paleontology, but it remains highly controversial.20,23,45 Various au- thors recognized the problems of inferring certain life history events, such as sexual maturity and primiparity, from growth functions.46 In taxo- nomic groups that do not have an extant analog, such as non-avian dinosaurs, arguments have been advanced in support of the use of the inflection point and of 90% asymptotic mass to infer ASM and age at onset of reproduction, respectively.46 However, this approach bears the intrinsic problem that, except for the Chapman-Richards model, the commonly used models (von Bertalanffy, Gompertz, and the logistic model) have a fixed inflection point, making the correct biological interpretation of this point problematic. Contrary to non-avian dinosaurs and other extinct phylogenetic groups, however, most mammals have extant analogs. This allows circumvention of the problematic use of a fixed inflection point. Instead, and as far as growth patterns and timing of life history events are known, these groups can be used as a phylo- genetically close extant taxon models for inferring ASM and other data in fossil species. iScience 26, 107654, September 15, 2023 5 ll iScience OPEN ACCESS Article Figure 5. Estimating ASM in N. rex (A) VB growth curves of L. europaeus (UnL): black dots; N. rex males: blue squares; N. rex females: red squares. SM: mean ASM of each sample. (B) Testing N. rex ASM deviation from allometric predictions using phylogenetic regression (PGLS). Residuals of ASM (y) against adult body mass (x) in log10 PGLS regression (Log10 Sexual Maturity = 1.762843 + 0,23227 * Log10 Body Mass). The points closest to N. rex are Marmota bobak (below N.rex to the right), insular Marmota vancouverensis (left to N. rex), and Erethizon dorsatum (below N.rex to the left; see also Figure S5). Hence, N. rex is comparable in size and age at maturity to these mammals; compared to extant leporids, however, it is a clear outlier (see discussion). Black dots: mammals (including ochotonids); orange dots: leporids; large green dot: M. balearicus; large red dot: N. rex. (Data from AnAge52). (C) Box-and-whisker plots depicting residuals from PGLS regression of ASM (days) against body mass (gr). Mammals (n = 84 taxa), leporids (n = 14 taxa), N. rex (n = 10), and M. balearicus. Horizontal lines denote the median, boxes the interquartile range (25–75 percent quartiles), and whiskers the maximum and minimum values; percentile ranks are computed by linear interpolation between the two nearest ranks. A one sample t test was performed (Past4) to compare the residuals of ASM with respect to BM mean of N. rex against the mammalian and leporid samples. The residual value of ASM (851 days) of N. rex differs significantly from the mammalian (t(82) = 35,074, 95% CI: 82.7–17.59, p = 6,3013E-50) and the leporid (t(14) = 46.77, 95% CI: 935.2–1025.1, p = 8,83E-17) means. The extant taxon model L. europaeus In the present study, our extant taxon model is based on the two subsamples: the FL subsample and the UnL subsample (see STAR Methods). In the FL subsample, age is controlled to the day by experimental protocols—accordingly, the FL growth model is based on accurate and detailed age data. However, the data available for N. rex consist of discrete years calculated from annual LAGs (skeletochronology); therefore, we needed to use the UnL subsample to obtain a growth model that is directly comparable with that of N. rex. We first checked the accuracy of the von Bertalanffy growth models (VBGM) obtained for the UnL and FL data. We calculated the statistical parameters of the two datasets, as well as those of the total Lepus sample (Tables S1-S3; Figures S1–S3). The results are similar in all instances (Table S4), validating the use of the UnL data to infer ASM in fossil leporids. Asymptotic size is a good estimator of ASM in leporids Using VBGM, we tested whether ASM in L. europaeus coincides with asymptotic size (i.e., 99% Linf; see STAR Methods), as suggested by empirical data.46 Determining the ASM of L. europaeus is not a simple question, considering that reproduction is interrupted in autumn and that it is not an instantaneous process. Broekhuizen & Maaskamp,47 using the presence of corpora lutea as the reproductive aptitude, indicate that most females become sexually mature at the age of 6–7 months, but they do not reproduce (i.e., get pregnant) until the age of 8–12 months. These results were later confirmed by the subsequent work of other authors.48–50 To calculate the ASM from VBGM, we used the method proposed by Ogle & Isermann51 to estimate age at a specific size from VBGM. UnL model provided ASM at 10.4 months, FL model at 8.5 months, and the total sample at 8.5 months, fitting empirical observations.46 For leporids, hence, age at 99% Linf is a good estimator of ASM (Figure 5A; Table S4). Inferring ASM in N. rex using L. europaeus as extant model Sexual dimorphism in N. rex We found two discrete groups of Nuralagus individuals with VBGMs (Tables S5-S7) that differ in asymptotic (Linf) values. The lack of overlap of the confidence intervals (Table S7) suggests weak but perceptible sexual dimorphism in the growth pattern of this taxon. 6 iScience 26, 107654, September 15, 2023 iScience ll Article OPEN ACCESS Figure 6. Epiphyseal fusion sequences in Lepus and Nuralagus Proximal and distal femora of N. rex (color) and L. europaeus (black and white); bottom left and right: slides from the respective femora of N. rex. Left: proximal and distal femora; from left to right: proximal posterior, proximal anterior, distal anterior, distal medial, distal lateral. First row: L. europaeus (IPS83881) 31 days; second row: N. rex right femur (IPS15586) R 3 years; third row: L. europaeus (IPS88711, male) 48 days; 4th row: L. europaeus (IPS84962) 79 days; bottom: slide of N. rex right femur (IPS15586) showing 3 winter LAGs (red) and a weaning LAG (green). Right: distal femora from left to right: anterior, medial, lateral; first row: L. europaeus (IPS 88858) 151 days; second row: N. rex right femur (IPS 11987, male) z 4 years; third row: L. europaeus (IPS 88859) 310 days; bottom: slide of N. rex right femur (IPS 11987, male; color inverted) showing 4 winter LAGs (red) and weaning LAG (green). Scales: 2 cm: all bones; 1000 mm: slides. The frequency distribution histogram of adults (mid-shaft transversal diameter of femora) shows two peaks that reflect these discrete groups (Figure S4), with the largest specimens having higher Linf than the smaller ones (see also Figure 5A). Females in sexual dimor- phic leporids are generally larger than males and have higher Linf.53 Attribution of larger N. rex specimens to females is further sup- ported by the fact that individuals with higher Linf stop growing earlier than those with lower Linf, the former being a typical trait of female mammals. ASM of N. rex To calculate ASM for N. rex, we applied the same model used for testing the UnL subsample of L. europaeus.51 The age at 99% of Linf was 3.6 years for larger and 6.2 years for smaller individuals (Figure 5A; Tables S5-S7). Taken together, the differences in asymptotic (Linf) values, the frequency distribution histogram of adults, and ASM suggest that the larger, earlier maturing N. rex individuals are females and the smaller, later maturing individuals are males. Following predictions from body mass scaling for herbivorous mammals (leporids, ochotonids, artiodactyls, and rodents up to 25 kg; AnAge database52), N. rex should have reached sexual maturity much earlier, at an age of 1.28 years (467 days). Considering that the N. rex sample has the same type of raw data as the UnL model (size: total surface area; age: annual lags), both samples are directly comparable. VBGM of N. rex strongly differs from that of L. europaeus (Figure 5A; Tables S3 and S4). Nuralagus’ growth constant k is extremely low (k = 0,63 for males; 1,08 for females) compared to that of L. europaeus (k = 5,06 for UnL subsample; k = 6,77 for FL subsample and k = 6,42 for the total sample). This indicates that N. rex growth rate was almost 5 times lower than that of L. europaeus (UnL sub-sample). Epiphyseal fusion sequences in Lepus and Nuralagus Epiphyseal fusion sequences are useful markers of maturity54; they allow identification and timing of equivalent maturational stages in closely related species (e.g.,55). Therefore, the state of epiphyseal fusion can be used to test the results of ASM from growth functions. Proximal femoral epiphyses fuse earlier than distal ones in both Nuralagus and Lepus Proximal epiphysis of the femoral head is already fusing at weaning in Lepus (IPS 83881, 31 days), whereas trochanter major remains unfused until approximately day 50 (IPS 88711, 48 days) (Figure 6). In Nuralagus, proximal femoral epiphysis starts fusing at the age of R3 years (after 3rd winter LAGs, Figure 6 bottom left, IPS 15586), long after weaning. At this age, trochanter major remains unfused. Thus, proximal epiph- yseal fusion in Nuralagus is importantly delayed and does not correspond to the same developmental stage as in Lepus. iScience 26, 107654, September 15, 2023 7 ll iScience OPEN ACCESS Article Table 1. Sample studied Species Sub-sample ID sex age origin Collection Lepus europaeus Unlabeled (UnL) IPS67842 f 5y 5m 2w FIWI ICP IPS67843 m 4y 5m FIWI ICP IPS67848 m 3y FIWI ICP IPS67858 f 3y FIWI ICP IPS67864 m 1y 1m FIWI ICP IPS67866 f 126d FIWI ICP IPS67870 m 240d FIWI ICP IPS67883 f 1y 3m FIWI ICP Labeled (FL) IPS82782 m 3d FIWI ICP IPS82783 f 11d FIWI ICP IPS82784 m 24d FIWI ICP IPS83878 f 8d FIWI ICP IPS83879 m 12d FIWI ICP IPS83880 f 28d FIWI ICP IPS83881 m 31d FIWI ICP IPS83882 m 34d FIWI ICP IPS84958 f 22d FIWI ICP IPS84959 m 54d FIWI ICP IPS84960 m 49d FIWI ICP IPS84961 m 57d FIWI ICP IPS84962 f 79d FIWI ICP IPS88707 m 3d FIWI ICP IPS88708 f 5d FIWI ICP IPS88709 m 12d FIWI ICP IPS88711 m 48d FIWI ICP IPS88712 m 102d FIWI ICP IPS88858 m 151d FIWI ICP IPS88859 m 309d FIWI ICP IPS93663 f 360d FIWI ICP Nuralagus rex IPS11987 m Punta Nati ICP IPS11989 m Punta Nati ICP IPS15163 f Punta Nati ICP IPS15375 juv R1.5y Punta Nati ICP IPS15518 f Punta Nati ICP IPS15586 f Punta Nati ICP IPS15629 f Punta Nati ICP IPS15679 m Punta Nati ICP IPS15334 f Punta Nati ICP IPS15574 f Punta Nati ICP IPS15951 f Punta Nati ICP IPS16100 f Punta Nati ICP ID: FIWI = Research Institute of Wildlife Ecology (Wien, Austria), ICP = Institut Català de Paleontologia Miquel Crusafont (Barcelona, Spain); sex: f/m: female/ male; age: d = days, w = weeks, m = months, y = years; origin: FIWI = Research Institute of Wildlife Ecology, Punta Nati (Menorca, Balearic Islands); collection: ICP = Institut Català de Paleontologia Miquel Crusafont (Barcelona, Spain). 8 iScience 26, 107654, September 15, 2023 iScience ll Article OPEN ACCESS Distal femoral epiphysis in Lepus starts fusing at an age between 50 and 80 days (Figure 6 IPS88711, 48 days; IPS84962, 79 days). By late puberty/beginning of sexual maturity (151 days/4.9 months), distal fusion is advanced and is completed shortly after beginning of sexual maturity46 (in our sample before 310 days/10.2 months, Figure 6 right, IPS88859). In Nuralagus, the suture of the distal epiphysis of the femur is still visible at the age of z4 years (IPS111987), which is comparable with the state of fusion in hares at beginning of sexual maturity (IPS88859). The age at distal epiphyseal fusion in Nuralagus is similar to the ASM esti- mated from growth curves as 3.6 and 6.2 years for females and males, respectively. Thus, N. rex maintains the leporid pattern of distal epiph- yseal fusion at an age around sexual maturity, supporting the results of ASM obtained from the growth model. Testing N. rex ASM deviation from allometric predictions using PGLS The island syndrome model suggests that long-term changes in insular traits, including delay of ASM, result from directional selection.1 For this notion to be true, ASM should be expected to (i) exceed the values predicted by allometric scaling and (ii) be independent of phylogeny. As we have seen before, bone histology, epiphyseal fusion, and VBGM provide evidence that giant Nuralagus attained sexual maturity far later than predicted from body mass scaling; it thus fulfills the first condition (i). We tested the second condition (ii) using phylogenetic generalized least squares regressions (PGLS), a proper framework to control for the potential influence of phylogeny (non-independence of the residuals).56 We furthermore implemented a comparative analysis of the residuals of ASM against body mass in a sample of 81 extant taxa of medium-sized (100–25.000 gr.) herbivorous mammals (Lagomorpha, Rodentia, Artiodactyla52; Figures S5; S6; Tables S8; S9; Data S1). Our results from phylogenetically corrected comparisons indicate that ASM in N. rex is clearly delayed beyond the values expected from body mass scaling in continental leporids and other mammals, and only slightly earlier than the extreme late ASM of the insular bovid M. balearicus17–20,57 (Figures 5 B, C; Figure S5; Table S8). Differences between residuals of N. rex and all continental taxa are statistically highly significant, particularly when compared with leporids (Figures 5 B, C). We thus conclude that the slow life history of the insular giant leporid does not merely result from allometric scaling but is an evolutionary response (adaptation) to the insular ecosystem as predicted by the island syndrome model and life history theory.1,4,7 DISCUSSION The importance of bone histology for understanding life history evolution on islands Within the last decades, bone histology has been shown to be a powerful tool to reliably reconstruct key life history traits in extant or subfossil insular vertebrates,16–21,23,25–29,57–59 where field data on the few surviving endemic vertebrate taxa are extremely difficult to obtain due to reduced population sizes, constraints imposed by conservation regulations for endangered species, and problems inherent to fieldwork. For this reason, the histological approach to the life history of fossil taxa that evolved unaffected by anthropogenic habitat destruction, intro- duced domestic animals, and other invasive mainland species is essential.17–21,23,25–29,57–59 It is precisely these endemics that provide deep insights into the undisturbed course of evolution on pristine islands. Though reproductive data cannot be reconstructed for N. rex, our histological comparisons with continental leporids provide evidence that the giant Menorcan leporid remarkably shifted toward the slow end of the reproductive-timing axis (sensu39). Indeed, residuals of ASM against body mass (Figure 5B) show that Nuralagus deviates from allometric predictions for extant leporids and is comparable only to a few similar-sized, extremely late-maturing mammals with decreased exposure to predators (low extrinsic mortality) such as large hiber- nating ground squirrels of the genus Marmota and Erethizon dorsatus, a spinous, arboreal hystricomorph from the North American continent. Marmota bobak, the steppe marmot, lives in the steppes of Eastern Europe and Central Asia, while M. vancouverensis, which has a similar age at maturity as N. rex, is an endemic insular taxon that lives in the high mountains of the Vancouver Island (Canada). The marked slowdown in growth rate, and significant delay in ASM compared to continental leporids and other continental mammalian herbivores of similar size ev- idence that Nuralagus is a clear outlier and, according to current knowledge, only surpassed by the dwarf insular bovid Myotragus,19,20,57 the only suitable fossil island mammal with available data within the taxonomic framework of our analysis. How do our results compare to other studies of insular giants? There are only very few studies that focus on age at maturity in insular giants.24–28,58–60 The results of all these studies suggest a delay in the age of maturity. However, none of these studies questioned whether the extension of the juvenile growth phase in these giants exceeded the predictions of allometric scaling models. One reason might be that an increase in age at maturity with body size is not surprising as it is pre- dicted by scaling rules.9,10 The other reason might be that the question has not aroused much interest since the growth pattern in small island vertebrates is not thought to be caused by resource scarcity as is the case in large insular vertebrates. The most common hypotheses to explain the evolution of gigantism in small island vertebrates are (i) decreased interspecific competition in depauperate insular ecosystems,61–65 (ii) release from predation pressure,66,67 (iii) high population densities and associated increased intra- specific competition,68 and (iv) local environmental conditions such as nutrient increase through guano and waste products of seabirds,68 ac- cess to aquatic resources,62 or elevated primary production sustained by high soil quality.69 As a result of one or more of these generally causally linked factors, insular small vertebrates are considered to have sufficient resources to evolve and maintain very large body sizes. Even more, they are deemed to attain their large size by maximizing growth rate.26,67 Seemingly iScience 26, 107654, September 15, 2023 9 ll iScience OPEN ACCESS Article paradoxically, however, a recent phylogenetic meta-analysis of paired island-mainland populations of terrestrial vertebrates found that small species are larger on islands with low productivity.67 While it is true that growing faster requires more resources, size increase can also be reached under resource constraints by growing over an extended time period. Indeed, the few studies that include both growth rate and growth period in their analysis have demonstrated that island giants grow longer, not faster, than their continental ancestors. Thus, Castanet and Baéz59 provided evidence that extant and subfossil Canarian giant lizards of the genus Gallotia do not accelerate growth; instead, they attain a larger size by growing at similar rates as smaller insular species, but over a longer period, thereby postponing the onset of reproduction. Hasegawa and Mori70 showed that young individuals of giant Japanese striped snake (Elaphe quadrivirgata) grow slower than conspecifics of other insular populations but achieve a larger asymp- totic size by doubling the growing time. The authors attribute the gigantic adult size to continued growth under low extrinsic mortality that results in increased longevity. Sandvig et al.,71 in a literature-based phylogenetic meta-analytic approach, compared growth rates in altricial insular, purportedly giant birds with those of mainland birds. They found a trend toward slower growth, which they suggest is associated with a delay in ASM, as predicted by Palkovacs’ theoretical life history model for age and size at maturity.7 The results of our PGLS analysis in an insular giant mammal provide the first evidence that insular gigantism is associated with a delay in age at maturity and a remarkably reduced growth rate k that far exceed predictions from allometry. The very slow and prolonged growth that eventually leads to an adult giant makes it clear that abundance of resources is not a prerequisite for gigantism. In agreement with life history theory, this solves the paradox67 that it is precisely on small, unproductive islands that small vertebrate taxa are evolving into giants. Our study underscores the importance of life history approaches and paleohistology as tools to elucidate fundamental questions about patterns and mechanisms in island evolution. Extensive empirical work, particularly in bone and dental histology, is still needed on extant and fossil island taxa to further test the generality of our results and to expand the still tenuous databases of important yet elusive life history traits of endangered long-lived island endemics. So far, our results confirm the life history axiom of the island syndrome model for giant insular mammals. LIMITATIONS OF THE STUDY Sample size is the most evident and pervasive limitation in paleontology. The number of specimens available is usually low. In histological studies that typically involve the destruction of the specimens, the number of available specimens is even more constrained. The quantity of the material of Nuralagus collected from the Menorcan site is limited, and the preparation of the bones—firmly embedded in the tough karstic sediment—was difficult and time-consuming. We made a huge effort to obtain the final number of 12 measurable individual femora. A further limitation is the narrow taxonomic focus. This has, at this stage, no solution because age and size at maturity, as well as the growth factor k (key life history characters here used to determine the position of island giants on the fast-slow life history axis), are not known for other fossil or extant insular giants. Given the enormous potential of hard-tissue histology and the inherent difficulties in collecting such data for largely protected island species, we propose to broaden the taxonomic focus and to include in future studies giants from other phylogenetic groups such as living cloud rats, the extinct Deinogalerix, the Sicilian giant dormouse Leithia, Caribbean Hutias, and giant sloths, to name just a few. The sample of extant European brown hares was constrained by time since it took us almost three years to label a useful number of in- dividuals of different ages and sexes to obtain a meaningful ontogenetic series. It is, therefore, an extraordinary sample that will certainly help solving important questions regarding the biology of L. europaeus as well as technical histological questions. STAR+METHODS Detailed methods are provided in the online version of this paper and include the following: d KEY RESOURCES TABLE d RESOURCE AVAILABILITY B Lead contact B Materials availability B Data and code availability d EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS B Ethical approval B Animal collection d METHODS DETAILS B Preparation of histological slides B Measurements B External fundamental system (EFS) d QUANTIFICATION AND STATISTICAL ANALYSIS B Von bertalanffy growth curve B Von bertalanffy growth curves using non-linear mixed effect models B Testing N. rex ASM deviation from allometric predictions using phylogenetic regression (PGLS) 10 iScience 26, 107654, September 15, 2023 iScience ll Article OPEN ACCESS SUPPLEMENTAL INFORMATION Supplemental information can be found online at https://doi.org/10.1016/j.isci.2023.107654. ACKNOWLEDGMENTS We thank Alessandro Urciuoli for his help with PGLS, M. Ángeles Meneses for her help with statistical inferences, and Manuel Fernandez for preparing the carcasses, and for making the histological sections of both the labeled and unlabeled specimens of European brown hares and the fossil specimens of N. rex. We would like to thank the three anonymous reviewers for their helpful and constructive comments. The study was supported by the R+D+I projects PID2020-117118GB-I00 (M.K., C.N.-M., J.Q.-C.) and PID2020-116908GB-I00 (S.M.-S.) funded by MCIN/ AEI/10.13039/501100011033/; by Generalitat de Catalunya (Centers de Recerca de Catalunya, CERCA Program); S.M.-S. received financial support by Generalitat de Catalunya (consolidated research groups SGR 86 and SGR 11686). The group ‘‘Paleoecology and Evolutionary Biology’’ is recognized without financial support by Generalitat de Catalunya (MK., CN-M SGR 01184). C.N.-M. is grant holder of the Programa Postdoctoral Beatriu de Pinós of the Secretaria d’Universitats i Recerca of the Generalitat de Catalunya (2021 00078 BP). The Research Institute of Wildlife Ecology provided additional financial and personnel support. AUTHOR CONTRIBUTIONS Conceptualization: M.K. Planning and Design: M.K. and S.M.S. Excavation and preparation of fossil material: J.Q.C. Implementation of the experimental part: W.A. and G.S; at an early stage also F.S. Methodology, Micrograph recordings: C.N.-M. and M.K.; Measurements: C.N.-M. and M.K.; Image editing: M.K; Histological analysis: M.K. and C.N.-M.; Statistical analysis: C.N.-M., S.M.-S. Writing: M.K. and S.M.-S.; C.N.-M. for parts of statistical descriptions. Discussion of the results: all authors. Funding Acquisition, M.K., S.M.-S., and W.A. DECLARATION OF INTERESTS The authors declare no competing interests. INCLUSION AND DIVERSITY We support inclusive, diverse, and equitable conduct of research. Received: April 7, 2023 Revised: May 23, 2023 Accepted: August 15, 2023 Published: August 17, 2023 REFERENCES 1. Adler, G.H., and Levins, R. (1994). The Island 11. Roth, V.L. (1990). Insular dwarf elephants: a 17. Köhler, M. (2009). The evolution of life history Syndrome in Rodent Populations. Q. Rev. case study in body mass estimation and traits associated to dwarfing in insular large Biol. 69, 473–490. ecological inference. In Body Size in mammals: a paleontological approach. 2. Foster, J.B. (1964). Evolution of mammals on Mammalian Paleobiology: Estimation and J. Vertebr. Paleontol. 29, 128A. islands. Nature 202. pp234–235. Biological Implications, J. Damuth and B. 18. Köhler, M. (2010). Fast or slow? The evolution 3. van Valen, L. (1973). A new evolutionary law. MacFadden, eds. (Cambridge University of life history traits associated with insular Evol. Theor. 1, 1–33. Press), pp. 151–180. dwarfing. In Islands and Evolution, V. Pérez- 4. Stearns, S.C. (1992). The Evolution of Life 12. Raia, P., Barbera, C., and Conte, M. (2003). Mellado and C. Ramon, eds. (Institut Histories (Oxford University Press). The fast life of a dwarfed giant. Evol. Ecol. 17, Menorquı́ d’Estudis. Recerca)), pp. 261–280. 5. Berrigan, D., and Koella, J.C. (1994). The 293–312. https://doi.org/10.1023/ 19. Köhler, M., and Moyà-Solà, S. (2009). evolution of reaction norms: simple models A:1025577414005. Physiological and life history strategies of a for age and size at maturity. J. Evol. Biol. 7, 13. Raia, P., and Meiri, S. (2006). The Island Rule in fossil large mammal in a resource-limited 549–566. https://doi.org/10.1046/j.1420- large mammals: Paleontology meets environment. Proc. Natl. Acad. Sci. USA. 106, 9101.1994.7050549.x. Ecology. Evolution 60, 1731–1742. https:// 20354–20358. https://doi.org/10.1073/pnas. 6. Stearns, S.C., and Koella, J.C. (1986). The doi.org/10.1554/05-664.1. 0813385106. evolution of phenotypic plasticity in life 14. Gould, S.J. (1977). Ontogeny and Phylogeny 20. Köhler, M., Herridge, V., Nacarino-Meneses, history traits: predictions of reaction norms (Belknap Press). C., Fortuny, J., Moncunill-Solé, B., Rosso, A., for age and size at maturity. Evol 40, 893–913. 15. Roth, V.L. (1992). Inferences from allometry Sanfilippo, R., Palombo, M.R., and https://doi.org/10.2307/2408752. and fossils: dwarfing of elephants on islands. Moyà-Solà, S. (2021). Palaeohistology reveals 7. Palkovacs, E.P. (2003). Explaining adaptive Oxf. Surv. Evol. Biol. 8, 259–288. a slow pace of life for the dwarfed Sicilian shifts in body size on islands: a life history 16. Bromage, T.G., Dirks, W., Erdjument- elephant. Sci. Rep. 11, 22862. approach. Oikos 103, 37–44. https://doi.org/ Bromage, H., Huck, M., Kullmer, O., Öner, R., 21. Kolb, C., Scheyer, T.M., Lister, A.M., Azorit, 10.1034/j.1600-0706.2003.12502.x. Sandrock, O., and Schrenk, F. (2002). A life C., de Vos, J., Schlingemann, M.A.J., Rössner, 8. Bonner, J.T. (1965). Size and Cycle: An Essay history and climate change solution to the G.E., Monaghan, N.T., and Sánchez-Villagra, on the Structure of Biology (Princeton evolution and extinction of insular dwarfs: a M.R. (2015). Growth in fossil and extant deer University Press). Cypriot experience. In World Islands in and implications for body size and life history 9. Peters, R.H. (1983). The Ecological Prehistory: International Insular evolution. BMC Evol. Biol. 15, 19. https://doi. Implications of Body Size (Cambridge Investigations, V Deia International org/10.1186/s12862-015-0295-3. University Press). Conference of Prehistory, W.H. Waldren and 22. Long, E.S., Courtney, K.L., Lippert, J.C., and 10. Calder, W.A.I. (1984). Size, Function and Life J.A. Ensenyat, eds. (Archaeopress)), Wall-Scheffler, C.M. (2019). Reduced body History (Dover Publications). pp. 420–427. size of insular black-tailed deer is caused by iScience 26, 107654, September 15, 2023 11 ll iScience OPEN ACCESS Article slowed development. Oecologia 189, 36. Palombo, M.R., Kohler, M., Moya Sola, S., and 739–757. https://doi.org/10.1007/s42991- 675–685. https://doi.org/10.1007/s00442- Giovinazzo, C. (2008). Brain versus body mass 021-00127-0. 019-04367-3. in endemic ruminant artiodactyls: A case 50. Stott, P., and Harris, S. (2006). Demographics 23. Hayashi, S., Kubo, M.O., Sánchez-Villagra, studied of Myotragus balearicus and smallest of the European hare (Lepus europaeus) in M.R., Izawa, M., Shiroma, T., Nakano, T., and Candiacervus species from Mediterranean the Mediterranean climate zone of Australia. Fujita, M. (2023). Variation and process of life Islands. Quat. Int. 182, 160–183. https://doi. Mamm. Biol. 71, 214–226. https://doi.org/10. history evolution in insular dwarfism as org/10.1016/j.quaint.2007.08.037. 1016/j.mambio.2006.02.009. revealed by a natural experiment. Front. 37. van der Geer, A.A.E., Lyras, G.A., Macphee, 51. Ogle, D.H., and Isermann, D.A. (2017). Earth Sci. 11, 1095903. https://doi.org/10. R.D.E., Lomolino, M., and Drinia, H. (2014). Estimating age at a specified length from the 3389/feart.2023.1095903. Mortality in a predator-free insular von Bertalanffy growth function. N. Am. J. 24. Todd, C.M., Westcott, D.A., Rose, K., Martin, environment: The dwarf deer of crete. Fish. Manag. 37, 1176–1180. https://doi.org/ J.M., and Welbergen, J.A. (2018). Slow Zootaxa 3847, 1–32. https://doi.org/10.1206/ 10.1080/02755947.2017.1342725. growth and delayed maturation in a Critically 3807.1. 52. Tacutu, R., Thornton, D., Johnson, E., Endangered insular flying fox (Pteropus 38. Hamada, F., and Mizuta, T. (2020). Unique Budovsky, A., Barardo, D., Craig, T., DIana, E., natalis). J. Mammal. 99, 1510–1521. https:// reproductive traits of the Amami rabbit Lehmann, G., Toren, D., Wang, J., et al. doi.org/10.1093/jmammal/gyy110. Pentalagus furnessi: an endangered endemic (2018). Human Ageing Genomic Resources: 25. Bourdon, E., Castanet, J., de Ricqlès, A., species from southwestern Japan. Mamm. New and updated databases. Nucleic Acids Scofield, P., Tennyson, A., Lamrous, H., and Res. 65, 805–813. https://doi.org/10.1007/ Res. 46, D1083–D1090. https://doi.org/10. Cubo, J. (2009). Bone growth marks reveal s13364-020-00497-9. 1093/nar/gkx1042. protracted growth in New Zealand kiwi (Aves, 39. Bielby, J., Mace, G.M., Bininda-Emonds, 53. Lindenfors, P., Gittleman, L.J., and Jones, E. Apterygidae). Biol. Lett. 5, 639–642. https:// O.R.P., Cardillo, M., Gittleman, J.L., Jones, (2007). Sexual size dimorphism in mammals. doi.org/10.1098/rsbl.2009.0310. K.E., Orme, C.D.L., and Purvis, A. (2007). The In Sex, Size and Gender Roles: Evolutionary 26. Angst, D., Chinsamy, A., Steel, L., and Hume, fast-slow continuum in mammalian life Studies of Sexual Size Dimorphism Sex, Size J.P. (2017). Bone histology sheds new light on history: An empirical reevaluation. Am. Nat. and Gender Roles: Evolutionary Studies of the ecology of the dodo (Raphus cucullatus, 169, 748–757. https://doi.org/10.1086/ Sexual Size Dimorphism, D.J. Fairbairn, W.U. Aves, Columbiformes). Sci. Rep. 7, 7993. 516847. Blanckenhorn, and S. Tamás, eds. (Oxford https://doi.org/10.1038/s41598-017-08536-3. 40. Kirkwood, T.B.L., and Shanley, D.P. (2005). University Press), pp. 1176–1180. 27. Chinsamy, A., Angst, D., Canoville, A., and Food restriction, evolution and ageing. 54. Brimacombe, C.S. (2017). The enigmatic Göhlich, U.B. (2020). Bone histology yields Mech. Ageing Dev. 126, 1011–1016. https:// relationship between epiphyseal fusion and insights into the biology of the extinct doi.org/10.1016/j.mad.2005.03.021. bone development in primates. Evol. elephant birds (Aepyornithidae) from 41. Jordana, X., Marı́n-Moratalla, N., Moncunill- Anthropol. 26, 325–335. https://doi.org/10. Madagascar. Biol. J. Linn. Soc. Lond. 130, Solè, B., Nacarino-Meneses, C., and Köhler, 1002/evan.21559. 268–295. https://doi.org/10.1093/biolinnean/ M. (2016). Ontogenetic changes in the 55. Brimacombe, C.S., Kuykendall, K.L., and blaa013. histological features of zonal bone tissue of Nystrom, P. (2018). Epiphyseal fusion and 28. de Ricqlès, A., Bourdon, E., Legendre, L.J., ruminants: a quantitative approach. C R dental development in Pan paniscus with and Cubo, J. (2016). Preliminary assessment Palevol 15, 255–266. https://doi.org/10.1016/ comparisons with Pan troglodytes. Am. J. of bone histology in the extinct elephant bird j.crpv.2015.03.008. Phys. Anthropol. 167, 903–913. https://doi. Aepyornis (Aves, Palaeognathae) from 42. Nacarino-Meneses, C., Jordana, X., and org/10.1002/ajpa.23710. Madagascar. C R Palevol 15, 197–208. https:// Köhler, M. (2016). Histological variability in 56. Smaers, J.B., and Rohlf, F.J. (2016). Testing doi.org/10.1016/j.crpv.2015.01.003. the limb bones of the Asiatic wild ass and its species’ deviation from allometric 29. Orlandi-Oliveras, G., Jordana, X., Moncunill- significance for life history inferences. PeerJ 4, predictions using the phylogenetic Solé, B., and Köhler, M. (2016). Bone e2580. https://doi.org/10.7717/peerj.2580. regression. Evolution 70, 1145–1149. https:// histology of the giant fossil dormouse 43. Nacarino-Meneses, C., and Chinsamy, A. doi.org/10.1111/evo.12910. Hypnomys onicensis (Gliridae, Rodentia) (2021). Mineralized-tissue histology reveals 57. Marı́n-Moratalla, N., Jordana, X., Garcı́a- from Balearic Islands. C R Palevol 15, protracted life history in the Pliocene three- Martı́nez, R., and Köhler, M. (2011). Tracing 238–244. https://doi.org/10.1016/j.crpv.2015. toed horse from Langebaanweg (South the evolution of fitness components in fossil 05.001. Africa). Zool. J. Linn. Soc. 196, 1117–1137, In bovids under different selective regimes. C R 30. Quintana, J., Köhler, M., and Moyà-Solà, S. press. https://doi.org/10.1093/zoolinnean/ Palevol 10, 469–478. https://doi.org/10.1016/ (2011). Nuralagus rex, gen. et sp. nov., an zlab037. j.crpv.2011.03.007. endemic insular giant rabbit from the 44. de Margerie, E., Cubo, J., and Castanet, J. 58. Guarino, F., and Andreone, F. (2003). Giant Neogene of Minorca (Balearic Islands, Spain). (2002). Bone typology and growth rate: and long-lived? Age structure in J. Vertebr. Paleontol. 31, 231–240. https:// testing and quantifying ‘Amprino’s rule’ in the Macroscincus coctei, an extinct skink from doi.org/10.1080/02724634.2011.550367. mallard (Anas platyrhynchos) C. C. R. Biol. Cape Verde. Amphib. Reptil. 24, 459–470. 31. Pons-Moyà, J., Moyà-Solà, S., Agustı́, J., and 325, 221–230. 59. Castanet, J., and Baéz, M. (1991). Adaptation Alcover, J.A. (1981). La fauna de mamı́feros 45. Werner, J., and Griebeler, E.M. (2014). and evolution in Gallotia lizards from the de los yacimientos menorquines con Allometries of maximum growth rate versus Canary Islands: age, growth, maturity and Geochelone symnesica (Bate, 1914). Nota body mass at maximum growth indicate that longevity. Amphib. Reptil. 12, 81–102. https:// preliminar. Acta geològica Hispànica 16, non-avian dinosaurs had growth rates typical doi.org/10.1163/156853891X00356. 129–130. of fast growing ectothermic sauropsids. PLoS 60. Turvey, S.T., Green, O.R., and Holdaway, R.N. 32. Moncunill-Solé, B., Jordana, X., Marı́n- One 9, e88834. https://doi.org/10.1371/ (2006). Cortical growth marks reveal Moratalla, N., Moyà-Solà, S., and Köhler, M. journal.pone.0088834. extended juvenile development in New (2014). How large are the extinct giant insular 46. Griebeler, E.M., Klein, N., and Sander, P.M. Zealand moa. Nature 435. https://doi.org/10. rodents? New body mass estimations from (2013). Aging, Maturation and Growth of 1038/nature03635. teeth and bones. Integr. Zool. 9, 197–212. Sauropodomorph Dinosaurs as Deduced 61. McNab, B.K. (1994). Resource use and the https://doi.org/10.1111/1749-4877.12063. from Growth Curves Using Long Bone survival of land and freshwater vertebrates on 33. Alcover, J., Moyà-Solà, S., and Pons-Moyà, J. Histological Data: An Assessment of oceanic islands. Am. Nat. 144, 643–660. (1981). Les quimeres del passat: els vertebrats Methodological Constraints and Solutions. 62. McNab, B.K. (2010). Geographic and fòssils del Plio-Quaternari de les Balears i PLoS One 8, 67012. https://doi.org/10.1371/ temporal correlations of mammalian size Pitiüses (Moll). journal.pone.0067012. reconsidered: A resource rule. Oecologia 34. Köhler, M., and Moyà-Solà, S. (2004). 47. Broekhuizen, S., and Maaskamp, F. (1981). 164, 13–23. https://doi.org/10.1007/s00442- Reduction of Brain and Sense Organs in the Annual production of young in European 010-1621-5. Fossil Insular Bovid Myotragus. Brain Behav. hares (Lepus europaeus) in the Netherlands. J 63. Meiri, S. (2007). Size evolution in island lizards. Evol. 63, 125–140. https://doi.org/10.1159/ Zool 193, 499–516. Glob. Ecol. Biogeogr. 16, 702–708. https:// 000076239. 48. Lelo, S., and Spahic, E. (2007). ZEC, Lepus doi.org/10.1111/j.1466-8238.2007.00327.x. 35. Köhler, M., and Moyà-Solà, S. (2001). europaeus Pallas, 1778 (Mammalia: 64. van der Geer, A., Lyras, G., de Vos, J., and Phalangeal adaptations in the fossil insular Lagomorpha). UZIZAZ 3, 46–51. Dermitzakis, M. (2010). Evolution of Island goat Myotragus. J. Vertebr. Paleontol. 21, 49. Heldstab, S.A. (2021). Habitat characteristics Mammals: Adaptation and Extinction of 621–624. https://doi.org/10.1671/0272- and life history explain reproductive Placental Mammals on Islands (Wiley- 4634(2001)021[0621:PAITFI]2.0.CO;2. seasonality in lagomorphs. Mamm. Biol. 101, Blackwell). 12 iScience 26, 107654, September 15, 2023 iScience ll Article OPEN ACCESS 65. Lomolino, M. v (1985). Body size of mammals 70. Hasegawa, M., and Mori, A. (2008). Does a 75. Rutzki, I., Elvers, W., Maisey, J., and Kellner, A. on islands: the Island Rule reexamined. Am. gigantic insular snake grow faster or live (1994). Chemical preparation techniques. In Nat. 125, 310–316. https://doi.org/10.1086/ longer to be gigantic? Evidence from a long- Vertebrate Paleontological Techniques, 284343. term field study. South Am. J. Herpetol. 3, Volume 1, P. Leiggi and P. May, eds. 66. Sondaar, P.Y. (1977). Insularity and its effect 145–154. https://doi.org/10.2994/1808- (Cambridge University Press), pp. 155–186. on mammal evolution. In Major patterns of 9798(2008)3[145:dagisg]2.0.co;2. 76. Nacarino-Meneses, C., and Köhler, M. (2018). vertebrate evolution, P.C. Goody, B.M. 71. Sandvig, E.M., Coulson, T., and Clegg, S.M. Limb bone histology records birth in Hecht, and M.K. Hecht, eds. (New York: (2019). The effect of insularity on avian growth mammals. PLoS One 13, e0198511. https:// Plenum), pp. 671–707. rates and implications for insular body size doi.org/10.1371/journal.pone.0198511. 67. Benı́tez-López, A., Santini, L., Gallego- evolution. Proc. Biol. Sci. 286, 20181967. Zamorano, J., Milá, B., Walkden, P., https://doi.org/10.1098/rspb.2018.1967. 77. Calderón, T., Arnold, W., Stalder, G., Painer, Huijbregts, M.A.J., and Tobias, J.A. (2021). 72. Orme, D., Freckleton, R., Thomas, G., J., and Köhler, M. (2021). Labelling The island rule explains consistent patterns of Petzoldt, T., Fritz, S.-S., Isaac, N., and experiments in red deer provide a general body size evolution in terrestrial vertebrates. Maintainer, W.P. (2018). Caper: Comparative model for early bone growth dynamics in Nat. Ecol. Evol. 5, 768–786. https://doi.org/ Analyses of Phylogenetics and Evolution in R. ruminants. Sci. Rep. 11, 14074. https://doi. 10.1038/s41559-021-01426-y. 73. Kawai, M., Delany, A.M., Green, C.B., Adamo, org/10.1038/s41598-021-93547-4. 68. Pafilis, P., Meiri, S., Foufopoulos, J., and M.L., and Rosen, C.J. (2010). Nocturnin 78. Lindstrom, M.L., and Bates, D.M. (1990). Valakos, E. (2009). Intraspecific competition Suppresses Igf1 Expression in Bone by Nonlinear Mixed Effects Models for and high food availability are associated with Targeting the 3 ’ Untranslated Region of Igf1 Repeated Measures Data. Biometrics 46, insular gigantism in a lizard. mRNA. Endocrinology 151, 4861–4870. 673–687. Naturwissenschaften 96, 1107–1113. https:// https://doi.org/10.1210/en.2010-0407. doi.org/10.1007/s00114-009-0564-3. 79. Bates, D.W., Mächler, M., Bolker, B.M., and 74. López-Martı́nez, N. (1989). Revisión Walker, S.C. (2015). Fitting linear mixed- 69. Meiri, S., Meijaard, E., Wich, S.A., Groves, sistemática y biostratigráfica de los effects models using lme4. BMJ Qual. Saf. 24, C.P., and Helgen, K.M. (2007). Mammals of Lagomorpha (Mammalia) del Terciario y 1–3. https://doi.org/10.18637/jss.v067.i01. Borneo - Small size on a large island. Cuaternario de España (Memorias del Museo J. Biogeogr. 35, 1087–1094. https://doi.org/ Paleontológico de Zaragoza. Diputación de 80. Wickham, H. (2016). ggplot2: Elegant 10.1111/j.1365-2699.2008.01897.x. Aragón)). Graphics for Data Analysis (Springer-Verlag). iScience 26, 107654, September 15, 2023 13 ll iScience OPEN ACCESS Article STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Experimental models: Organisms/strains Lepus europaeus See Table 1 See Table 1 Fossils of Nuralagus rex See Table 1 See Table 1 Chemicals, peptides and recombinant proteins Alizarin complexone VWR International GmbH Cat#1.01010.0005 Calcein green VWR International GmbH Cat#1.02315.0005 Ketamine Richter Pharma AG, Austria Ketamidor 100 mg/ml Medetomidine Le Vet B.V, the Netherlands Narcostart 1mg/ml T 61 injection solution for euthanasia Intervet, Vienna, Austria Embutramide, Mebezonium and Tetracaine Software and algorithms Lme4 package in R (v. 4.1.0) R Core Team, 2020 http://www.r-project.org/ Phylogenetic generalized least-square regression (PGLS) R software Orme et al.72 RESOURCE AVAILABILITY Lead contact Further information and requests should be directed to the lead contact, Meike Köhler (meike.kohler@icp.cat). Materials availability Histological slides generated in this study belong to the collections of the Institut Català de Palaeontologia and can be studied upon request. Data and code availability d All data reported in this paper will be shared by the lead contact upon request. d This paper does not report original code. d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS Ethical approval All experiments and procedures were approved by the Ethics and Animal Welfare Committee of the University of Veterinary Medicine, Vienna in accordance with the University’s guidelines for Good Scientific Practice and authorized by the Austrian Federal Ministry of Education, Sci- ence and Research (permit number: BMWFW-68.205/0101-WF/II/3b/2014) in accordance with current legislation. Furthermore, these exper- iments were carried out in compliance with ARRIVE guidelines. Animal collection L. europaeus We used a large sample of L. europaeus as an extant comparison taxon because: I. hares nicely show the formation of all tissue types and developmental stages within a short time span as they grow very rapidly (as we observed in a previous incomplete labeling experiment with hares from the same population at the Research Institute of Wildlife Ecol- ogy (FIWI) II. Weaning in hares occurs comparably rapidly as several litters can be born in the same year. This makes it possible to focus on the most interesting growth periods and examine multiple individuals in detail. We subdivided the hares into two groups: (i) A fluorochrome-labelled subsample (FL) and (ii) an un-labelled subsample (UnL). 14 iScience 26, 107654, September 15, 2023 iScience ll Article OPEN ACCESS Fluorochrome-labelled subsample (FL). All markers were administered subcutaneously. Alizarin was administered at a dose of 30 mg/kg, and Calcein green at a dose of 8 mg/kg. Solutions were prepared under sterile conditions at the pharmacy of the Vetmeduni Vienna and buffered to 7.4 pH with NaHCO3. We labeled the significant number of 24 hares of both sexes (2 does, 16 leverets, (+6 back up leverets), of which data from 21 in- dividuals were usable (Table 1). Back up animal were only enrolled if needed and with an empirically estimated general mortality of less than thirty per cent. The FL individuals were labelled from birth (some as early as 3 days before birth and still intrauterine) to about one year of age (past ASM). Leverets were weaned early (day 27–28) or late (day 35–44) and are from either spring, summer or fall litters in order to detect possible differences in bone accretion rate between them. After tetracycline was discarded because antibiotics should no longer be used in labelling experiments, and Alizarin red and Xylenol Orange were excluded due to their toxicity, we started using the following dyes: Calcein Green, Alizarin complexone, and Calcein blue. However, as we observed under the microscope, cal- cein blue faded almost instantly. Thus, we only used Calcein Green and Alizarin complexone. We always injected the markers at the same time of day (2:00 p.m.) to avoid potential changes in bone apposition rate that might be induced by circadian Igf1 transcript levels.73 Animals were sacrificed at almost any age to preserve bone tissue that would otherwise be lost at later ages through enlarge- ment of the medullary cavity. We established a high-resolution labelling schedule to generate a reliable sequence of ontogenetic stages up to ASM. This allowed cal- ibrating early life history events (birth and weaning) with bone histology and detecting potential differences in growth rate and expression of tissue growth markers generated by differences in season of birth and timing of weaning. Animals were managed by keepers with decades of European brown hare experience. Veterinary care was constantly available and guar- anteed by veterinarians working at the Research Institute of Wildlife Ecology (FIWI), who have extensive experience in management of hares. Does were kept in the breeding population of the FIWI. Animals are usually group housed in co-sexual groups of up to four animals in 6x6 m floor enclosures enriched with straw bedded hiding boxes, branches and pawing boxes. After mating, does were transferred to indi- vidual standard breeding boxes (2x1 m) that are structured with 50% solid flooring with straw/nesting material and hiding space and 50% pervious, open flooring. The does could be separated from their offspring via grafters, to allow them to retreat after suckling, similar to natural conditions. All enclosures are in the FIWI stable in which animals are exposed to natural light but protected from direct sunlight. Food and water were available to all animals ad libitum at all times. Animals were checked daily, cages were cleaned and supplied with fresh bedding two times a week. Body mass was measured daily. From about 120 days onwards, BM was measured over periods of 5 days, each followed by 2-3 days of interruption. For euthanasia, animals were deeply anesthetized with a combination of Ketamine (30 mg/kg,) and Medetomidine (250 mg/kg), in succes- sion they received an IC overdose of Embutramide, Mebezonium and Tetracaine (T 61, Intervet, Vienna Austria). Unlabelled subsample (UnL). The second subsample of unlabelled individuals consists of 8 specimens of known sex and age at death and was used for skeletochronology, based on annual winter LAGs (Table 1). This sample extends the timeframe beyond ASM to 5 years of age. Individuals were also kept at FIWI during their life, in similar conditions as those described for the FL subsample. The fossil menorcan leporid N. rex We used 12 femoral slides of N. rex to calculate the growth curves from the ontogenetic series of the successive annual (winter) LAGs. Of these, 8 were identified as females and 3 as males; and a juvenile specimen (Table 1). Preparation of nuralagus fossils. The remains of the vertebrates of the Punta Nati karstic deposits are included in very hard red silt strongly cemented with calcium carbonate (Figure 1C). The removal and recovery of the bones from the matrix was carried out in a aqueous solution of acetic acid at a concentration of 10 % following the criteria of López-Martı́nez74 and Rutzky et al.,.75 Before the first immersion in the acid solution, the bones located on the rock surface were protected and consolidated with acryloid B67. Once dried, the samples were subjected to acid attack for two days. After that time, the blocks were placed in buckets with freshwater for two days in order to eliminate the remnants of acetic acid and calcium acetate. The remains included in the red clay sediments depos- ited at the bottom of the buckets were washed on a very fine sieve. The obtained concentrate was also immersed in freshwater for two days. Once the water bath was finished, the samples were dried, and the bones were consolidated again. In order to streamline the process, two groups of samples were worked on simultaneously: while one was treated with acid, the other one was dried and consolidated. The acid treat- ment lasted until the rock blocks were completely dissolved, which took place for around four years. The bone concentrates, recovered by sieving, were kept apart until their final processing (physical separation and identification). Once separated from the matrix, the surfaces of the larger bones had numerous impurities due to a mixture of acryloid and red silt. These residues were carefully removed, taking into account the fragility of the bones due to their decalcification. The process ended once the bones were completely clean and properly consolidated. iScience 26, 107654, September 15, 2023 15 ll iScience OPEN ACCESS Article METHODS DETAILS Preparation of histological slides After euthanasia, the carcasses of hares were frozen in the FIWI’s refrigeration facilities and flown to the Institut Català de Paleontologia (ICP) by a specialized transport company. At the ICP, all skeletal parts were prepared using cold water baths (so as not to destroy the label colours) with a weakly concentrated KH-7 degreasing agent and stored in the ICP collections. Histological slides were prepared according to our laboratory’s standard protocol.29,76,77 The long bones were cut and a chunk of about 2 cm was extracted from the mid-diaphysis of each bone. The blocks were degreased and dehydrated by alcohol immersion and embedded in Araldite 2020 epoxy resin. Then they were cut in half with a low-speed diamond saw (IsoMet low speed saw, Buehler). The cut surfaces were polished with a MetaServ polishing machine and fixed with epoxy resin on a frosted glass. Once the sample was fixed, it was cut to a thickness of 100–120 mm with a diamond saw (Petrothin, Buehler) and finely polished again to obtain the finished slide. Histological thin sections were observed under a Zeiss Scope.A1 microscope, labelled samples with fluorescence filters. The recordings were made with the built-in camera (AxioCam ICc5). Measurements The data for the von Bertalanffy growth curves come from measurements on micrographs of bone cross-sections using a Zeiss Scope.A1 mi- croscope with optional fluorescence and polarization settings, and an integrated camera (AxioCam ICc5). For Nuralagus and Lepus (UnL) we traced each annual LAG and measured the entire cross-sectional area bordered by the LAG including the medullary cavity with Photoshop (2023, 24.4.1). We included the medullary cavity in the total area to avoid additional noise in the dataset from age- and diet-related variations in the size and shape of the medullary cavity resulting from increased bone resorption and the deposition of new bone in the inner cortical surface in response to changing biomechanical stresses and strains over the course of ontogeny. This way, the increasing values of consec- utive LAGs represent the respective total area for each year. Each bone series thus corresponds to the ontogenetic growth trajectory of an individual. For labeled Lepus (FL), the labels were traced and, the same as previously described for LAGs, their area including the medullary cavity was measured, giving day-to-day size/age increments. Again, each bone series corresponds to an individual We propose replacing the commonly used radial distances 44 with area measurements as here described. With a constant regular increase in bone surface area, radial distances become progressively shorter, visually simulating a decrease in cortical growth rate (which can be cor- rected by dividing the distance by the time elapsed between the two limiting LAGs 44). Furthermore, local differences in tissue growth due to biomechanical stresses and strains can lead to locally very different radial distances, potentially affecting the modelling of growth curves. Increasing expansion of transverse bone surface area is not subject to local fluctuations; transverse bone surface area at midshaft tightly cor- relates body weight and can be taken as its proxy for this. Because of these considerations, the use of areas instead of lengths is especially appropriate in studies of slow-growing mammals and those with irregular bone cortical outline. External fundamental system (EFS) Identifying an EFS is not difficult in hares (see Figures 2A and 2B). Its starts around the end of puberty with the onset of reproductive maturity, which we found in the femora of both individuals, the female (IPS93663a-1 ID16) and the late-born male (IPS88859 ID15) to be shortly after 8 months, beginning with the first winter LAG (beginning of February - beginning of March). Thus, the pubertal lamellar tissue with simple vascularization that starts around 5 months of age, does not form part of the EFS, although in histological studies it is usually included in this tissue structure because of the similarity with the EFS tissue. In hares, hence, the EFS begins with the first winter LAG. The EFS is a rather thin layer even in old animals, because radial growth resumes only weakly in the following spring and stops in almost all individuals after 2-3 years. Because of the extremely slow growth, tissue deposition might cease several times during winter leading to several thin LAGs that appear as one clear winter LAG each year. The ESF is the last growth stage before complete arrest (no further LAGs). The ESF in Nuralagus is much more complicated to identify (see Figure 4; female IPS15574). In fact, the femoral tissue of Nuralagus grows extremely slowly. In the first winter, there is a transition from lamellar tissue with sparse primary osteons to dense lamellar tissue with infre- quent, scattered simple vascular canals (Figures 3 and 4). In the following (2nd) year, the tissue is almost avascular, and tissues of subsequent years are completely avascular. Thus, the year-to-year tissue transitions are subtle and almost imperceptible because of the extremely slow growth rate, which lets the overall tissue from the first winter onwards appear pubertal / mature. Nevertheless, the radial distances between the first three LAGs indicate that this slow-down is not comparable with the (quasi) cessation of tissue apposition in mature hares (Figure 4). Indeed, it is not before the fourth winter that annual LAGs are hard to count, suggesting that puberty may have started just before, or around the third winter of life when the vascularization disappeared completely, and reproductive maturity began between the third and the fourth winter (Figure 4). This coincides with the results from von Bertalanffy growth curves based on bone cross sectional areas. QUANTIFICATION AND STATISTICAL ANALYSIS Von bertalanffy growth curve The Von Bertalanffy growth curve is a widely used growth model in which individual growth is treated as an increase in length or weight with increasing age (in this study cross sectional bone area / year). The typical von Bertalanffy growth curve is represented by the equation  y = Linf  1  eð  Kðt  t0Þ 16 iScience 26, 107654, September 15, 2023 iScience ll Article OPEN ACCESS where: y is the length or weight (in our case bone area, see ’Measurements’), t is time (in our case years) and Linf, K, and t0 are the growth param- eters. Linf represents the asymptotic length or weight (in our case bone area). K is also called the Brody growth rate coefficient, although it is not a real growth rate. It measures the exponential rate of approach to the asymptotic size. Finally, t0 is a model artifact the adjust the model for the initial size of the animal.51 Von bertalanffy growth curves using non-linear mixed effect models Nonlinear mixed-effects models are subject-specific models in which a general growth model is assumed to characterize the population, but the coefficients of the growth model can be unique to the individual. This statistical model incorporates both fixed effects and random effects. Fixed effects are population parameters assumed to be the same each time data is collected, and random effects are random variables asso- ciated with each sample (individual) from a population. Longitudinal data, such as those in our study (i.e., more than two LAGs or chemical labels per individual), can be modeled using non-linear mixed-effect models.78,79 We used lme4 package79 in R software (v. 4.1.0) to adjust our different datasets (Nuralagus Large, Nuralagus Small, Lepus All, Lepus LAGs and Lepus Labelled) to a mixed-effect von Bertalanffy model via the maximum likelihood method. Akaike’s Information Criterion (AIC), Bayesian Information Criterion (BIC), the log-likelihood value (LogLik) and the coefficient of determination (R2) were used to assess the ac- curacy of the fitted models. For each dataset, we estimate an initial fit in which we added random effects on the three parameters of the model (Linf, K, t0). If the variance of the random effect for any of the parameters was below 0.001, we performed a second fit in which we eliminate the random effect for this specific parameter. We also allow for possible correlation of the random effects. The different fitted models were compared using ANOVA. We used the von Bertalanffy equations obtained for each dataset to estimate the end of growth, i.e., the age at maturity. Specifically, we resolved the different equations for y=99%Linf to obtain the time for this life history trait. All graphs were created using ggplot2 package80 in R software (v. 4.1.0). Testing N. rex ASM deviation from allometric predictions using phylogenetic regression (PGLS) Comparative allometric analysis of ASM and body mass in medium size (100 -25.000 gr.) herbivorous mammals (Lagomorphs, Rodentia, Ar- tiodactyla) of a sample of 81 extant taxa (AnAge database,52; Table S9), provides the best framework to assess if the N. rex ASM corresponds to expectations from body size scaling or there is a delay in maturity beyond expectations from scaling, as Palkovacs’7 theoretical model pre- dicts. To test the impact of body mass on the age at sexual maturity, while controlling for potential influence of phylogeny (non-independence of the residuals), the phylogenetic generalized least-square (PGLS) regression was performed.72 iScience 26, 107654, September 15, 2023 17