Philippa Kaur

Monitoring sites and speciesThe monitoring survey was conducted at five sites in Fukushima Prefecture (sites 1–4 and A1) and at one site in Ibaraki Prefecture (site 5), Japan (Fig. 1). Sites 1, 2, and A1 are located in Kawauchi Village, site 3 in Otama Village, site 4 in Tadami Town, and site 5 in Ishioka City. Monitoring at sites 1–5 was started in 2011 or 2012, and site A1 was additionally monitored since 2017. The tree species, age, mean diameter at breast height, initial deposition density of 137Cs, and sampling year of each sample at each site are listed in Table 1. The dominant tree species in the contaminated area, namely, Japanese cedar (Cryptomeria japonica [L.f.] D.Don), Japanese cypress (Chamaecyparis obtusa [Siebold et Zucc.] Endl.), konara oak (Quercus serrata Murray), and Japanese red pine (Pinus densiflora Siebold et Zucc.) were selected for monitoring. Japanese chestnut (Castanea crenata Siebold et Zucc.) was supplementally added in 2017. The cedar, cypress, and pine are evergreen coniferous species, and the oak and chestnut are deciduous broad-leaved species. Sites 1 and 3 each have three plots, and each plot contains a different monitoring species. Site A1 has one plot containing two different monitoring species, and the remaining sites each have one plot with one monitoring species, giving ten plots in total.Figure 1Locations of the monitoring sites and initial deposition densities of 137Cs (decay-corrected to July 2, 2011) following the Fukushima nuclear accident in Fukushima and Ibaraki Prefectures. Open circles indicate the monitoring sites and the cross mark indicates the Fukushima Dai-ichi Nuclear Power Plant. Data on the deposition density were provided by MEXT19,20 and refined by Kato et al.21. The map was created using R (version 4.1.0)22 with ggplot2 (version 3.3.5)23 and sf (version 1.0–0)24 packages.Full size imageTable 1 Description of the sampled trees and monitoring sites.Full size tableSample collection and preparationBulk sampling of bark and wood disks was conducted by felling three trees per year at all sites during 2011–20168,25 and at sites 3–5 and A1 during 2017–2020. Partial sampling from six trees per year was conducted at sites 1 and 2 during 2017–2020 (from seven trees at site 2 in 2017) to sustain the monitoring trees. All the samples were obtained from the stems around breast height. During the partial sampling, bark pieces sized approximately 3 cm × 3 cm (axial length × tangential length) were collected from four directions of the tree stem using a chisel, and 12-mm-diameter wood cores were collected from two directions of the tree stem using an automatic increment borer (Smartborer, Seiwa Works, Tsukuba, Japan) equipped with a borer bit (10–101-1046, Haglöf Sweden, Långsele, Sweden). Such partial sampling increases the observational errors in the bark and wood 137Cs concentrations in individual trees26. To mitigate this error and maintain an accurate mean value of the 137Cs concentration, we increased the number of sampled trees from three to six. The sampling was conducted mainly in July–September of each year; the exceptions were site-5 samples in 2011 and 2012, which were collected irregularly during January–February of the following year. The collected bark pieces were separated into outer and inner barks, and the wood disks and cores were split into sapwood and heartwood. The outer and inner bark samples during 2012–2016 were obtained by partial sampling of barks sized approximately 10 cm × 10 cm from 2–3 directions on 2–3 trees per year.The bulk samples of bark, sapwood, and heartwood were air-dried and then chipped into flakes using a cutting mill with a 6-mm mesh sieve (UPC-140, HORAI, Higashiosaka, Japan). The pieces of the outer and inner bark were chipped into approximately 5 mm × 5 mm pieces using pruning shears, and the cores of the sapwood and heartwood were chipped into semicircles of thickness 1–2 mm. Each sample was packed into a container for radioactivity measurements and its mass was measured after oven-drying at 75 °C for at least 48 h. Multiplying this mass by the conversion factor (0.98 for bark and 0.99 for wood)8 yielded the dry mass at 105 °C.Radioactivity measurementsThe radioactivity of 137Cs in the samples was determined by γ-ray spectrometry with a high-purity Ge semiconductor detector (GEM20, GEM40, or GWL-120, ORTEC, Oak Ridge, TN). For measurements, the bulk and partial samples were placed into Marinelli containers (2.0 L or 0.7 L) and cylindrical containers (100 mL or 5 mL), respectively. The peak efficiencies of the Marinelli containers, the 100-mL container, and the 5-mL container were calibrated using standard sources of MX033MR, MX033U8PP (Japan Radioisotope Association, Tokyo, Japan), and EG-ML (Eckert & Ziegler Isotope Products, Valencia, CA), respectively. For the measurement of the 5-mL container, a well-type Ge detector (GWL-120) was used under the empirical assumption that the difference in γ-ray self-absorption between the standard source and the samples is negligible27. The measurement was continued until the counting error became less than 5% (higher counting errors were allowed for small or weakly radioactive samples). The activity concentration of 137Cs in the bark (whole) collected by partial sampling was calculated as the mass-weighted mean of the concentrations in the outer and inner barks; meanwhile, the concentration in the wood (whole) was calculated as the cross-sectional-area-weighted mean of sapwood and heartwood concentrations. The activity concentrations were decay-corrected to September 1, 2020, to exclude the decrease due to the radioactive decay.Trend analysesThe yearly representative values (true states) of 137Cs activity concentration in each stem part in each plot were estimated using a DLM, a state-space model in which the noise follows a normal distribution and the relationship between variables is linear. One basic DLM is the local linear trend model defined by the following equations:$$Y_{t} = mu _{t} + varepsilon _{t} ,quad quad quad varepsilon _{t} sim Normal left( {0,sigma _{varepsilon }^{2} } right)$$
(1)
$$mu_{t} = mu_{t – 1} + beta_{t – 1} + eta_{t} ,quad quad quad eta_{t} sim Normal left( {0,sigma_{eta }^{2} } right)$$
(2)
$$beta_{t} = beta_{t – 1} + zeta_{t} ,quad quad quad zeta_{t} sim Normal left( {0,sigma_{zeta }^{2} } right)$$
(3)
where Yt, μt, and βt are the observation values, level (true state), and slope, respectively, and εt, ηt, and ζt denote their corresponding noises. The subscript t is the time index. The noises εt, ηt, and ζt follow normal distributions with a mean of 0 and variances of ({sigma }_{varepsilon }^{2}), ({sigma }_{eta }^{2}), and ({sigma }_{zeta }^{2}), respectively. To detect relatively long-term trends, we employed the smooth local linear trend model28 (also called the smooth trend model, integrated random walk model, or second-order trend model), which is obtained by considering that μt and βt are driven by the same noise. The trend changes are assumed to be smoother in this model than in the local linear trend model28,29. Combining Eqs. (2) and (3), μt in the smooth local linear trend model is finally obtained as$$mu_{t} = 2mu_{t – 1} – mu_{t – 2} + eta_{t} ,quad quad quad eta_{t} sim Normal left( {0,sigma_{eta }^{2} } right)$$
(4)
The parameters μt, ({sigma }_{eta }^{2}), and ({sigma }_{varepsilon }^{2}) of each stem part in each plot were determined by Bayesian estimation with a Markov chain Monte Carlo (MCMC) method. The Bayesian estimation was performed in R (version 4.1.0)22 with the rstan package (version 2.21.2)30. Uninformative prior distributions were used for μ1, μ2, ({sigma }_{eta }^{2}), and ({sigma }_{varepsilon }^{2}). The log-transformed values of the 137Cs activity concentration (decay-corrected to September 1, 2020) were given as Yt (the observed values of multiple individuals in each year were passed via the segment function of Stan). MCMC sampling was conducted for four chains of 50,000 iterations (the first 25,000 were discarded as warmup), obtaining 100,000 MCMC samples for each parameter. The MCMC was judged to have converged when the maximum value of Rhat was less than 1.05 and the divergent transitions after warmup were fewer than 1,000 (i.e., less than 1% of the MCMC sample size). On the datasets of the outer and inner barks from site-3 oaks and all stem parts from site-A1 pines and chestnuts, the MCMC converged poorly owing to the small number of monitoring years. Thus, the temporal trends in these datasets were not analyzed (the observational data at site A1 are shown in Supplementary Fig. S1 and Table S1).To detect decadal trends rather than yearly variations, we determined the temporal trends in the true state (μ) by setting 2–4 delimiting years and examining whether μ varied significantly from one delimiting year to the next. As the delimiting years, we selected the initial and final years of monitoring and the years in which the median µ was highest (µ-max year) and lowest (µ-min year). When the µ-max year and/or the µ-min year coincided with the initial year and/or final year of monitoring, the number of delimiting years reduced from four to two or three. The trend in µ between two delimiting years was determined to be increasing and decreasing when the 95% credible interval of µ2nd delimiting year − µ1st delimiting year (obtained from the MCMC samples) was higher and lower than zero, respectively. A flat trend (no significant variation) was detected when the 95% credible interval included zero. If the 3rd and 4th delimiting years existed, the trends between the 2nd and 3rd delimiting years and between the 3rd and 4th delimiting years were determined in the same manner.The 137Cs CRs of outer bark/inner bark, heartwood/sapwood, and inner bark/sapwood were also subjected to the above trend analyses. On datasets with less than five years of monitoring, the MCMC did not converge so the trend analysis was not attempted. More

Wild animals construct various types of structures that are adaptive to their life and reproduction. For example, termites that inhabit the African savanna use soil to construct a huge mound that reaches 10 m in height; they produce hollows and holes in these mounds to allow air ventilation, thereby keeping the internal temperature constant1. In addition, prairie dogs inhabiting the North American prairie dig vertically and horizontally extending burrows in the ground that they use for shelter and rearing offspring; these burrows have multiple entrances, some of which are chimney-shaped to improve ventilation efficiency2. In the field of biomimetics, researchers apply the principles of animal-created structures in applications useful to humans3.The white-spotted pufferfish Torquigener albomaculosus (Pisces: Tetraodontidae) is a relatively small species that grows to ~10 cm in total total length (Fig. 1). Male T. albomaculosus individuals construct an intricate geometric circular structure, known as the “mystery circle,” with a diameter of 2 m in the sand of the seabed;4 the discovery of these structures has fascinated researchers and the general public worldwide. The male pufferfish digs the sand on the seabed with its fins and body while swimming straight ahead toward the centre from different directions, and a circular structure composed of radially aligned peaks and valleys was constructed. Finally, the male creates a maze-like pattern by flapping its anal fin on the bottom of the central zone4. Thus, the male completes the circular structure by himself. Furthermore, we discovered that the earliest stage of the mystery circle is composed of dozens of irregular depressions, which might function as landmarks for the formation of the radial patterns5. By accumulating observations of pufferfish behaviour, we were able to conduct a computer simulation including the swimming trajectory of the pufferfish extracted from video images wherein they constructed the circular structure. This simulation revealed that an elaborate circular geometric pattern is inevitably formed if the pufferfish repeats the digging behavior on the seabed using simple rules6. We also observed the reproductive behaviour of the pufferfish and found that they consistently breed in a semilunar cycle from spring to summer. Each male constructs a mystery circle and spawns with multiple females on the nest, and the male cares for the eggs alone until they hatch. Some of the elements of the circular structure, i.e., its size, symmetry, ornaments, and maze-like pattern, might be important factors in terms of female mate choice4,7.Fig. 1The white-spotted pufferfish Torquigener albomaculosus. Lateral view of a male (a), and male digging behaviour on the seabed while rolling up fine sand particles (b).Full size imageAlthough data on the reproductive ecology and circle-construction behaviour of these pufferfish have been collected, many questions remain. Our interdisciplinary research currently has two themes: (i) theoretical studies on the logic of 3D-structure formation of the circular structure and (ii) ethological studies on the relationship between female mate choice and the features of the structure. To advance these studies, it is essential to collect quantitative data on the circular structure. Thus, we reconstructed 3D models of six completed mystery circles using a “structure from motion” (SfM) algorithm (Fig. 2).Fig. 2“Mystery circle” constructed by a white-spotted pufferfish (Torquigener albomaculosus). 3D model displayed on a computer (a), one of the video frames used to reconstruct the 3D model (b), and a Styrofoam model output in full size created using a 3D printer and the 3D data (c) for a specific mystery circle 20160615_K13.Full size imageOn the other hand, the mystery circle constructed by the pufferfish may have potential applications in biomimetics similar to the structures constructed by termites and prairie dogs. To support the importance of its structural characteristics, it has been observed that the water passing through the valley upstream always gathers in the center of the structure, regardless of the direction of water flow4. Furthermore, particle size analysis of the sand forming the mystery circle has revealed that it has the function of extracting fine-grained sand particles from the valleys arranged radially to the outside and directing them to the center (Kawase, in prep.). The field of computational fluid dynamics, which makes full use of fluid dynamics technology, engineering knowledge, and computers, will logically clarify the characteristics of the 3D structure of the mystery circle we have reconstructed here. Shameem et al. reconstructed a 3D model of a mystery circle to explore the flow features with 2D computational fluid dynamic simulations8. Since our model has already been quantified as 3D data, computational fluid analysis can be immediately performed using this data, and the structural features of the mystery circle are expected to be applied in a wide range of fields, such as architecture and engineering, via biomimetics. More

Butchart, S. H. M. et al. Global biodiversity: Indicators of recent declines. Science 328(5982), 1164–1168. https://doi.org/10.1126/science.1187512 (2010).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Dirzo, R. et al. Defaunation in the anthropocene. Science 345(6195), 401–406. https://doi.org/10.1126/science.1251817 (2014).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Bradshaw, C. J. A. et al. Underestimating the challenges of avoiding a ghastly future. Front. Conserv. Sci. https://doi.org/10.3389/fcosc.2020.615419 (2021).Article 

Google Scholar 
Barnosky, A. D. et al. Has the Earth’s sixth mass extinction already arrived?. Nature 471(7336), 51–57. https://doi.org/10.1038/nature09678 (2011).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Ceballos, G. et al. Accelerated modern human-induced species losses: Entering the sixth mass extinction. Sci. Adv. https://doi.org/10.1126/sciadv.1400253 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Ceballos, G., Ehrlich, P. R. & Raven, P. H. Vertebrates on the brink as indicators of biological annihilation and the sixth mass extinction. Proc. Natl. Acad. Sci. U.S.A. 117(24), 13596–13602. https://doi.org/10.1073/pnas.1922686117 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
McGowan, P. J., Traylor-Holzer, K. & Leus, K. IUCN guidelines for determining when and how ex situ management should be used in species conservation. Conserv. Lett. 10(3), 361–366. https://doi.org/10.1111/conl.12285 (2016).Article 

Google Scholar 
Clout, M. N. & Merton, D. V. Saving the Kakapo: The conservation of the world’s most peculiar parrot. Bird Conserv. Int. 8(3), 281–296. https://doi.org/10.1017/s0959270900001933 (1998).Article 

Google Scholar 
Milinkovitch, M. C. et al. Genetic analysis of a successful repatriation programme: Giant Galápagos tortoises. Proc. R. Soc. B Biol. Sci. 271(1537), 341–345. https://doi.org/10.1098/rspb.2003.2607 (2004).CAS 
Article 

Google Scholar 
Ryder, O. A. & Wedemeyer, E. A. A cooperative breeding programme for the Mongolian wild horse Equus przewalskii in the United States. Biol. Conserv. 22(4), 259–271. https://doi.org/10.1016/0006-3207(82)90021-0 (1982).Article 

Google Scholar 
Mallinson, J. J. C. Conservation breeding programmes: An important ingredient for species survival. Biodivers. Conserv. 4(6), 617–635. https://doi.org/10.1007/bf00222518 (1995).Article 

Google Scholar 
Seddon, P. J., Armstrong, D. P. & Maloney, R. F. Developing the science of reintroduction biology. Conserv. Biol. 21(2), 303–312. https://doi.org/10.1111/j.1523-1739.2006.00627.x (2007).Article 
PubMed 

Google Scholar 
Bowkett, A. E. Recent captive-breeding proposals and the return of the ark concept to global species conservation. Conserv. Biol. 23(3), 773–776. https://doi.org/10.1111/j.1523-1739.2008.01157.x (2009).Article 
PubMed 

Google Scholar 
Shan, L. et al. Large-scale genetic survey provides insights into the captive management and reintroduction of giant pandas. Mol. Biol. Evol. 31(10), 2663–2671. https://doi.org/10.1093/molbev/msu210 (2014).CAS 
Article 
PubMed 

Google Scholar 
Fischer, J. & Lindenmayer, D. An assessment of the published results of animal relocations. Biol. Conserv. 96(1), 1–11. https://doi.org/10.1016/s0006-3207(00)00048-3 (2014).Article 

Google Scholar 
Christie, M. R., Marine, M. L., French, R. A. & Blouin, M. S. Genetic adaptation to captivity can occur in a single generation. Proc. Natl. Acad. Sci. U.S.A. 109(1), 238–242. https://doi.org/10.1073/pnas.1111073109 (2011).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Fraser, D. J. et al. Population correlates of rapid captive-induced maladaptation in a wild fish. Evol. Appl. 12(7), 1305–1317. https://doi.org/10.1111/eva.12649 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Ralls, K., Brugger, K. & Ballou, J. Inbreeding and juvenile mortality in small populations of ungulates. Science 206(4422), 1101–1103. https://doi.org/10.1126/science.493997 (1979).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Charlesworth, D. & Charlesworth, B. Inbreeding depression and its evolutionary consequences. Annu. Rev. Ecol. Evol. Syst. 18(1), 237–268. https://doi.org/10.1146/annurev.es.18.110187.001321 (1987).Article 

Google Scholar 
Ralls, K., Ballou, J. D. & Templeton, A. Estimates of lethal equivalents and the cost of inbreeding in mammals. Conserv. Biol. 2(2), 185–193. https://doi.org/10.1111/j.1523-1739.1988.tb00169.x (1988).Article 

Google Scholar 
Hedrick, P. W. & Kalinowski, S. T. Inbreeding depression in conservation biology. Annu. Rev. Ecol. Evol. Syst. 31(1), 139–162. https://doi.org/10.1146/annurev.ecolsys.31.1.139 (2000).Article 

Google Scholar 
Frankham, R. Introduction to Conservation Genetics 2nd edn. (Cambridge University Press, 2010).Book 

Google Scholar 
Laikre, L. Conservation genetics of Nordic carnivores: Lessons from zoos. Hereditas 130(3), 203–216. https://doi.org/10.1111/j.1601-5223.1999.00203.x (2004).Article 

Google Scholar 
Gomendio, M., Cassinello, J. & Roldan, E. R. S. A comparative study of ejaculate traits in three endangered ungulates with different levels of inbreeding: Fluctuating asymmetry as an indicator of reproductive and genetic stress. Proc. R. Soc. B Biol. Sci. 267(1446), 875–882. https://doi.org/10.1098/rspb.2000.1084 (2000).CAS 
Article 

Google Scholar 
Swinnerton, K. J., Groombridge, J. J., Jones, C. G., Burn, R. W. & Mungroo, Y. Inbreeding depression and founder diversity among captive and free-living populations of the endangered pink pigeon Columba mayeri. Anim. Conserv. 7(4), 353–364. https://doi.org/10.1017/s1367943004001556 (2004).Article 

Google Scholar 
Farquharson, K. A., Hogg, C. J. & Grueber, C. E. Offspring survival changes over generations of captive breeding. Nat. Commun. https://doi.org/10.1038/s41467-021-22631-0 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Kleiman, D. G., Thompson, K. V. & Baer, C. K. Wild Mammals in Captivity: Principles and Techniques for Zoo Management 2nd edn. (University of Chicago Press, 2021).
Google Scholar 
Ralls, K. & Ballou, J. D. Captive breeding and reintroduction. In Encyclopedia of Biodiversity (ed. Levin, S. A.) 662–667 (Academic Press, 2013). https://doi.org/10.1016/b978-0-12-384719-5.00268-9.Chapter 

Google Scholar 
Reed, D. H. & Frankham, R. Correlation between fitness and genetic diversity. Conserv. Biol. 17(1), 230–237. https://doi.org/10.1046/j.1523-1739.2003.01236.x (2003).Article 

Google Scholar 
Spielman, D., Brook, B. W. & Frankham, R. Most species are not driven to extinction before genetic factors impact them. Proc. Natl. Acad. Sci. U.S.A. 101(42), 15261–15264. https://doi.org/10.1073/pnas.0403809101 (2003).ADS 
CAS 
Article 

Google Scholar 
Willi, Y., van Buskirk, J. & Hoffmann, A. A. Limits to the adaptive potential of small populations. Annu. Rev. Ecol. Evol. Syst. 37(1), 433–458. https://doi.org/10.1146/annurev.ecolsys.37.091305.110145 (2006).Article 

Google Scholar 
Habel, J. C., Husemann, M., Finger, A., Danley, P. D. & Zachos, F. E. The relevance of time series in molecular ecology and conservation biology. Biol. Rev. 89(2), 484–492. https://doi.org/10.1111/brv.12068 (2013).Article 
PubMed 

Google Scholar 
Araki, H., Cooper, B. & Blouin, M. S. Genetic effects of captive breeding cause a rapid, cumulative fitness decline in the wild. Science 318(5847), 100–103. https://doi.org/10.1126/science.1145621 (2007).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Purohit, D. et al. Genetic effects of long-term captive breeding on the endangered pygmy hog. PeerJ 9, e12212. https://doi.org/10.7717/peerj.12212 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Hahn, E. E. & Culver, M. Genetic diversity and structure in Arizona pronghorn following conservation efforts. Conserv. Sci. Pract. https://doi.org/10.1111/csp2.498 (2021).Article 

Google Scholar 
Charlesworth, B. Effective population size and patterns of molecular evolution and variation. Nat. Rev. Gen. 10(3), 195–205. https://doi.org/10.1038/nrg2526 (2009).CAS 
Article 

Google Scholar 
Wang, J., Santiago, E. & Caballero, A. Prediction and estimation of effective population size. Heredity 117(4), 193–206. https://doi.org/10.1038/hdy.2016.43 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
O’Gara, W., Yoakum, J. D. & McCabe, R. E. Pronghorn: Ecology and Managment (University Press of Colorado, 2004).
Google Scholar 
Janis, C. M., Scott, K. M. & Jacobs, L. L. Evolution of Tertiary Mammals of North America: Terrestrial Carnivores, Ungulates, and Ungulate like Mammals Vol. 1 (Cambridge University Press, 2005).
Google Scholar 
Nelson, E. W. Status of the Pronghorn Antelope, 1922–1924 (U.S Department Agriculture Bulletin, 1925).Book 

Google Scholar 
O’Gara, B. W. & McCabe, R. E. From exploitation to conservation. In Pronghorn: Ecology and Management (eds O’Gara, B. W. & Yoakum, J. D.) 41–73 (University Press Colorado, 2004).
Google Scholar 
Cancino, J., Ortega-Rubio, A. & Sanchez-Pacheco, J. A. Status of an endangered subspecies: The peninsular pronghorn at Baja California. J. Arid Environ. 32(4), 463–467. https://doi.org/10.1006/jare.1996.0039 (1996).ADS 
Article 

Google Scholar 
Laliberte, A. S. & Ripple, W. J. Range contractions of North American carnivores and ungulates. Bioscience 54(2), 123–138. https://doi.org/10.1641/0006-3568 (2004).Article 

Google Scholar 
Medellín, R. A. et al. History, ecology, and conservation of the pronghorn antelope, bighorn sheep, and black bear in Mexico. In Biodiversity, Ecosystems, and Conservation in Northern Mexico (eds Cartron, J.-L. et al.) 387–405 (Oxford University Press, 2005).
Google Scholar 
Lee, T. E., Bickham, J. W. & Scott, M. D. Mitochondrial DNA and allozyme analysis of North American pronghorn populations. J. Wildl. Manag. 58(2), 307–318. https://doi.org/10.2307/3809396 (1994).Article 

Google Scholar 
IUCN SSC Antelope Specialist Group. Antilocapra americana ssp. peninsularis. The IUCN Red List of Threatened Species 2021: e.T1679A200726719. https://doi.org/10.2305/IUCN.UK.2021-2.RLTS.T1679A200726719.en (2021).SEMARNAT. Norma Oficial Mexicana NOM-059-SEMARNAT-2010, Protección ambiental– Especies nativas de México de flora y fauna silvestres– Categorías de riesgo y especificaciones para su inclusión, exclusión o cambio– Lista de especies en riesgo. Diario Oficial de la Federación 30 diciembre (2010).U. S. Fish and Wildlife Service. Recovery Plan for the Sonoran pronghorn (Antilocapra americana sonoriensis), Second Revision. (U.S. Fish and Wildlife Service, Southwest Region, Albuquerque, 2016).Cancino, J., Sanchez-Sotomayor, V. & Castellanos, R. From the field: Capture, hand-raising, and captive management of peninsular pronghorn. Wildl. Soc. Bull. 33(1), 61–65. https://doi.org/10.2193/0091-7648 (2005).Article 

Google Scholar 
Horne, J. S., Hervert, J. J., Woodruff, S. P. & Mills, L. S. Evaluating the benefit of captive breeding and reintroductions to endangered Sonoran pronghorn. Biol. Conserv. 196, 133–146. https://doi.org/10.1016/j.biocon.2016.02.005 (2016).Article 

Google Scholar 
CONANP. Programa de Acción para la Conservación de la Especie: Berrendo (Antilocapra americana), 2009 año del berrendo. Secretaria del Medio Ambiente y Recursos Naturales (SEMARNAT). www.conanp.gob.mx (2009).Cancino, J., Rodríguez-Estrella, R. & Miller, P. Using population viability analysis for management recommendations of the endangered endemic peninsular pronghorn. Acta Zool. Mex. 26(1), 173–189 (2010).
Google Scholar 
Danoff-Burg, J. A. & Mulroe, K. Peninsular Pronghorn Species Action Plan (2021) (in press).Stephen, C. L. et al. Population genetic analysis of sonoran pronghorn (Antilocapra americana sonoriensis). J. Mammal. 86(4), 782–792. https://doi.org/10.1644/1545-1542 (2005).Article 

Google Scholar 
Stephen, C. L., Whittaker, D. G., Gillis, D., Cox, L. L. & Rhodes, O. E. Genetic consequences of reintroductions: An example from oregon pronghorn antelope (Antilocapra americana). J. Wildl. Manag. 69(4), 1463–1474. https://doi.org/10.2193/0022-541x (2005).Article 

Google Scholar 
Barnow-Meyer, K. & Byers, J. Genetic diversity and gene flow in Yellowstone Basin pronghorn (Antilocapra americana). UW Natl. Parks Serv. Res. Station Annu. Rep. 31, 65–72. https://doi.org/10.13001/uwnpsrc.2008.3705 (2008).Article 

Google Scholar 
LaCava, M. E. F. et al. Pronghorn population genomics show connectivity in the core of their range. J. Mammal. 101(4), 1061–1071. https://doi.org/10.1093/jmammal/gyaa054 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Klimova, A., Munguia-Vega, A., Hoffman, J. I. & Culver, M. Genetic diversity and demography of two endangered captive pronghorn subspecies from the Sonoran Desert. J. Mammal. 95(6), 1263–1277. https://doi.org/10.1644/13-mamm-a-321 (2014).Article 

Google Scholar 
Hahn, E. E., Klimova, A., Munguía-Vega, A., Clark, K. B. & Culver, M. Use of museum specimens to refine historical pronghorn subspecies boundaries. J. Wildl. Manag. 84(3), 524–533. https://doi.org/10.1002/jwmg.21810 (2020).Article 

Google Scholar 
Axelrod, D. I. The evolution of desert vegetation in western North America. Carnegie Instit. Wash. Publ. 590, 215–306 (1950).
Google Scholar 
Dolby, G. A., Bennett, S. E. K., Lira-Noriega, A., Wilder, B. T. & Munguía-Vega, A. Assessing the geological and climatic forcing of biodiversity and evolution Surrounding the Gulf of California. J. Southwest. 57, 391–455. https://doi.org/10.1353/jsw.2015.0005 (2015).Article 

Google Scholar 
Gedir, J. V., Cain, J. W., Harris, G. & Turnbull, T. T. Effects of climate change on long-term population growth of pronghorn in an arid environment. Ecosphere 6(10), art189. https://doi.org/10.1890/es15-00266.1 (2015).Article 

Google Scholar 
Cornuet, J. M. et al. DIYABC v2.0: A software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics 30(8), 1187–1189. https://doi.org/10.1093/bioinformatics/btt763 (2014).CAS 
Article 
PubMed 

Google Scholar 
Islas-Espinoza, M. & de las Heras, A. Peninsular pronghorn conservation: Too many paradigms, too few indicators. In Sustainability Indicators in Practice (eds Latawiec, A. & Agol, D.) 126–145 (De Gruyter Open Poland, 2015). https://doi.org/10.1515/9783110450507-012.Chapter 

Google Scholar 
Willoughby, J. R. et al. The impacts of inbreeding, drift and selection on genetic diversity in captive breeding populations. Mol. Ecol. 24(1), 98–110. https://doi.org/10.1111/mec.13020 (2014).CAS 
Article 
PubMed 

Google Scholar 
Crow, J. F. & Kimura, M. An Introduction in Population Genetics Theory (Harper and Row, 1970).MATH 

Google Scholar 
Falconer, D. S. Introduction to Quantitative Genetics 3rd edn. (Longman Scientific and Technical, 1989).
Google Scholar 
Ballou, J. D. Strategies for maintaining genetic diversity in captive populations through reproductive technology. Zoo Biol. 3(4), 311–323. https://doi.org/10.1002/zoo.1430030404 (1984).Article 

Google Scholar 
Ballou, J. D. & Lacy, R. C. Identifying genetically important individuals for management of genetic diversity in pedigreed populations. In Population Management for Survival and Recovery (eds Ballou, J. D. et al.) 76–111 (Columbia Press, 1995).
Google Scholar 
Montgomery, M. E. et al. Minimizing kinship in captive breeding programs. Zoo Biol. 16(5), 377–389. https://doi.org/10.1002/(sici)1098-2361 (1997).Article 

Google Scholar 
Dunn, S. J., Clancey, E., Waits, L. P. & Byers, J. A. Inbreeding depression in pronghorn (Antilocapra americana) fawns. Mol. Ecol. 20(23), 4889–4898. https://doi.org/10.1111/j.1365-294x.2011.05327.x (2011).Article 
PubMed 

Google Scholar 
Hoffman, J. I. et al. High-throughput sequencing reveals inbreeding depression in a natural population. Proc. Natl. Acad. Sci. U.S.A. 111(10), 3775–3780. https://doi.org/10.1073/pnas.1318945111 (2014).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Kardos, M. et al. The crucial role of genome-wide genetic variation in conservation. Proc. Natl. Acad. Sci. U.S.A. 118(48), e2104642118. https://doi.org/10.1073/pnas.2104642118 (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Zoonomia Consortium. A comparative genomics multitool for scientific discovery and conservation. Nature 587(7833), 240–245. https://doi.org/10.1038/s41586-020-2876-6 (2020).ADS 
CAS 
Article 

Google Scholar 
Ceballos, F. C., Joshi, P. K., Clark, D. W., Ramsay, M. & Wilson, J. F. Runs of homozygosity: Windows into population history and trait architecture. Nat. Rev. Genet. 19(4), 220–234. https://doi.org/10.1038/nrg.2017.109 (2018).CAS 
Article 
PubMed 

Google Scholar 
Supple, M. A. & Shapiro, B. Conservation of biodiversity in the genomics era. Genome Biol. https://doi.org/10.1186/s13059-018-1520-3 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Hohenlohe, P. A. & Rajora, O. P. Population Genomics: Wildlife (Springer, 2020).
Google Scholar 
Chalmers, G. A. & Barrett, M. W. Capture myopathy in pronghorns in Alberta, Canada. J. Am. Vet. Med. Assoc. 171(9), 918–923 (1977).CAS 
PubMed 

Google Scholar 
Sotelo-Gallardo, H., Contreras Balderas, A. J. & Espinosa Treviño, A. Comparación de dos métodos de liberación del berrendo, Antilocapra americana (Artiodactyla: Antilocapridae) en Coahuila, México. Rev. Biol. Trop. 65(3), 1208. https://doi.org/10.15517/rbt.v65i3.29447 (2017).Article 

Google Scholar 
Breed, D. et al. Conserving wildlife in a changing world: Understanding capture myopathy—A malignant outcome of stress during capture and translocation. Conserv. Physiol. https://doi.org/10.1093/conphys/coz027 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Snyder, N. F. et al. Limitations of captive breeding in endangered species recovery. Conserv. Biol. 10(2), 338–348. https://doi.org/10.1046/j.1523-1739.1996.10020338.x (1996).Article 

Google Scholar 
Bonebrake, T. C., Christensen, J., Boggs, C. L. & Ehrlich, P. R. Population decline assessment, historical baselines, and conservation. Conserv. Lett. 3(6), 371–378. https://doi.org/10.1111/j.1755-263x.2010.00139.x (2010).Article 

Google Scholar 
Grismer, L. L. & McGuire, J. A. The oases of central Baja California, Mexico. Part I. A preliminary account of the relict mesophilic herpetofauna and the status of the oases. Bull. South. Calif. Acad. Sci. 92, 2–24 (1993).
Google Scholar 
Welsh, H. H., Clark, W. H., Franco-Vizcaíno, E. & Valdéz-Villavicencio, J. H. Herpetofauna associated with palm oases across the Californian-Sonoran transition in Northern Baja California, Mexico. Southwest. Nat. 55(4), 581–585. https://doi.org/10.1894/pas-15.1 (2010).Article 

Google Scholar 
Mann, D. H., Groves, P., Gaglioti, B. V. & Shapiro, B. A. Climate-driven ecological stability as a globally shared cause of Late Quaternary megafaunal extinctions: The Plaids and Stripes Hypothesis. Biol. Rev. 94(1), 328–352. https://doi.org/10.1111/brv.12456 (2018).Article 

Google Scholar 
Brown, D. E., Warnecke, D. & McKinney, T. Effects of midsummer drought on mortality of doe pronghorn (Antilocapra americana). Southwest. Nat. 51(2), 220–225. https://doi.org/10.1894/0038-4909 (2006).Article 

Google Scholar 
Simpson, D. C., Harveson, L. A., Brewer, C. E., Walser, R. E. & Sides, A. R. Influence of precipitation on pronghorn demography in Texas. J. Wildl. Manag. 71(3), 906–910. https://doi.org/10.2193/2005-753 (2007).Article 

Google Scholar 
McKinney, T., Brown, D. E. & Allison, L. Winter precipitation and recruitment of pronghorns in Arizona. Southwest. Nat. 53(3), 319–325. https://doi.org/10.1894/cj-147.1 (2008).Article 

Google Scholar 
Otte, A. Partners save the Sonoran pronghorn. Endang. Species Bull. 31, 22–23 (2006).
Google Scholar 
McCullough, D. R. & Barrett, R. H. Wildlife 2001: Populations (Springer, 1992).Book 

Google Scholar 
Percie Du Sert, N. et al. Reporting animal research: Explanation and elaboration for the ARRIVE guidelines 20. PLOS Biol. 18(7), e3000411. https://doi.org/10.1371/journal.pbio.3000411 (2020).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Carling, M. D., Passavant, C. W. & Byers, J. A. DNA microsatellites of pronghorn (Antilocapra americana). Mol. Ecol. Not. 3(1), 10–11. https://doi.org/10.1046/j.1471-8286.2003.00334.x (2002).Article 

Google Scholar 
Dunn, S. J. et al. Ten polymorphic microsatellite markers for pronghorn (Antilocapra americana). Conserv. Genet. Resour. 2(1), 81–84. https://doi.org/10.1007/s12686-009-9166-9 (2010).Article 

Google Scholar 
Munguia-Vega, A., Klimova, A. & Culver, M. New microsatellite loci isolated via next-generation sequencing for two endangered pronghorn from the Sonoran Desert. Conserv. Genet. Resour. 5(1), 125–127. https://doi.org/10.1007/s12686-012-9749-8 (2012).Article 

Google Scholar 
Boutin-Ganache, I., Raposo, M., Raymond, M. & Deschepper, C. F. M13-Tailed primers improve the readability and usability of microsatellite analyses performed with two different allele-sizing methods. Biotechniques 31(1), 25–28. https://doi.org/10.2144/01311bm02 (2001).Article 

Google Scholar 
Amos, W. et al. Automated binning of microsatellite alleles: Problems and solutions. Mol. Ecol. Not. 7(1), 10–14. https://doi.org/10.1111/j.1471-8286.2006.01560.x (2006).CAS 
Article 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2021). https://www.R-project.org/.Jombart, T. & Ahmed, I. Adegenet 1.3–1: New tools for the analysis of genome-wide SNP data. Bioinformatics 27(21), 3070–3071. https://doi.org/10.1093/bioinformatics/btr521 (2011).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Kamvar, Z. N., Tabima, J. F. & Grünwald, N. J. Poppr: An R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ 2, e281. https://doi.org/10.7717/peerj.281 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
Adamack, A. T. & Gruber, B. PopGenReport: Simplifying basic population genetic analyses in R. Methods Ecol. Evol. 5(4), 384–387. https://doi.org/10.1111/2041-210x.12158 (2014).Article 

Google Scholar 
Agapow, P. M. & Burt, A. Indices of multilocus linkage disequilibrium. Mol. Ecol. Not. 1(1–2), 101–102. https://doi.org/10.1046/j.1471-8278.2000.00014.x (2001).CAS 
Article 

Google Scholar 
Paradis, E. pegas: An R package for population genetics with an integrated-modular approach. Bioinformatics 26(3), 419–420. https://doi.org/10.1093/bioinformatics/btp696 (2010).CAS 
Article 
PubMed 

Google Scholar 
Goudet, J. hierfstat, a package for r to compute and test hierarchical F-statistics. Mol. Ecol. Not. 5(1), 184–186. https://doi.org/10.1111/j.1471-8286.2004.00828.x (2005).Article 

Google Scholar 
Aparicio, J. M., Ortego, J. & Cordero, P. J. What should we weigh to estimate heterozygosity, alleles or loci?. Mol. Ecol. 15(14), 4659–4665. https://doi.org/10.1111/j.1365-294x.2006.03111.x (2006).CAS 
Article 
PubMed 

Google Scholar 
Alho, J. S., Välimäki, K. & Merilä, J. Rhh: An R extension for estimating multilocus heterozygosity and heterozygosity–heterozygosity correlation. Mol. Ecol. Res. 10(4), 720–722. https://doi.org/10.1111/j.1755-0998.2010.02830.x (2010).Article 

Google Scholar 
Stoffel, M. A. et al. inbreedR: An R package for the analysis of inbreeding based on genetic markers. Methods Ecol. Evol. 7(11), 1331–1339. https://doi.org/10.1111/2041-210x.12588 (2016).Article 

Google Scholar 
Wang, J. Coancestry: A program for simulating, estimating and analysing relatedness and inbreeding coefficients. Mol. Ecol. Res. 11(1), 141–145. https://doi.org/10.1111/j.1755-0998.2010.02885.x (2010).ADS 
Article 

Google Scholar 
Wang, J. Triadic IBD coefficients and applications to estimating pairwise relatedness. Genet. Res. 89(3), 135–153. https://doi.org/10.1017/s0016672307008798 (2007).CAS 
Article 
PubMed 

Google Scholar 
Marshall, T. C. et al. Estimating the prevalence of inbreeding from incomplete pedigrees. Proc. R. Soc. B Biol. Sci. 269(1500), 1533–1539. https://doi.org/10.1098/rspb.2002.2035 (2002).CAS 
Article 

Google Scholar 
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. https://doi.org/10.18637/jss.v067.i01 (2015).Article 

Google Scholar 
Beaumont, M. A., Zhang, W. & Balding, D. J. Approximate Bayesian computation in population genetics. Genetics 162(4), 2025–2035. https://doi.org/10.1093/genetics/162.4.2025 (2002).Article 
PubMed 
PubMed Central 

Google Scholar 
Bertorelle, G., Benazzo, A. & Mona, S. ABC as a flexible framework to estimate demography over space and time: Some cons, many pros. Mol. Ecol. 19(13), 2609–2625. https://doi.org/10.1111/j.1365-294x.2010.04690.x (2010).CAS 
Article 
PubMed 

Google Scholar 
Fagundes, N. J. R. et al. Statistical evaluation of alternative models of human evolution. Proc. Natl. Acad. Sci. U.S.A. 104(45), 17614–17619. https://doi.org/10.1073/pnas.0708280104 (2007).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar  More

Sample collection and hydrographySampling was conducted aboard the University of Victoria’s MSV John Strickland either weekly, biweekly or monthly between 11 March 2010 and 15 November 2017 in Saanich Inlet at 48.59°N, 123.50°W (Fig. 1). To standardize measurements and due to biological significance, seawater was collected from the euphotic zone. Sampling depths corresponded to approximately 100, 50, 15, and 1% of the photosynthetically active radiation (PAR) at the surface (Io). These “light” depths were either determined using a CTD-mounted PAR sensor or a Secchi disk. CTD profiles were performed prior to each seawater cast to measure depth, temperature and conductivity of the water column, and PAR, fluorescence, and dissolved oxygen (when available).Seawater from each light depth was collected using Niskin or GO-FLO bottles on either a rosette sampler or an oceanographic wire. When possible, individual samples were collected directly from the Niskin or GO-FLO bottles. When time was not sufficient to allow direct sampling, bulk samples of seawater from each depth were collected into acid-washed polyethylene carboys, kept cold in the dark, and homogenized before sub-sampling for the individual measurements.Dissolved nutrientsFor the measurements of nitrite (NO2−), nitrate and nitrite (NO3− + NO2−), phosphate (PO43−) and Si(OH)4, seawater samples from each light depth were syringe-filtered through a combusted 0.7 µm (nominal porosity) glass fibre filter into acid-washed 30-mL polypropylene bottles and immediately frozen. All nutrient samples were stored at −20 °C until analysis. Concentrations of NO2−, NO3− + NO2−, PO43−, and Si(OH)4 were determined using an Astoria Nutrient Autoanalyzer (Astoria-Pacific, OR, USA) following the methodology of Barwell-Clarke and Whitney22. During 2014 and 2015, samples for the measurement of Si(OH)4 were collected separately from those for the other nutrients, filtered with a 0.6 µm polycarbonate membrane filter and stored at 4 °C. During this period, Si(OH)4 concentrations were determined manually using the molybdate blue colorimetric methodology23. Replicate (2 or 3) nutrient samples were taken at each depth; average data are presented in the published dataset and the figures (Fig. 2).Fig. 2Dissolved macronutrient concentrations in the euphotic zone of Saanich Inlet from March 11, 2010 to November 15, 2017. Left panels show depth-integrated concentrations (black bars on top) and time-series profiles (filled contour/scatter plots on bottom) for (A) nitrate plus nitrite (NO3− + NO2−), (B) phosphate (PO43−) and (C) silicic acid (Si(OH)4). In the time-series profiles, 2012–2013 data are not interpolated due to single-depth sampling. Grey shaded regions in top panels indicate the phytoplankton growing seasons considered for this study (March 1st – October 30th). Right panels show monthly-averaged depth profiles for the entire 8-year period, illustrating euphotic zone seasonality for each nutrient. The color scale bars on the far right apply to both the time-series vertical profiles and the 8-year seasonal plots. Sampling depths are indicated by round symbols. The year labels are positioned under the tick marks corresponding to January.Full size imageSuspended particulate matterTotal chlorophyll-aChlorophyll-a (Chl-a) was used as a proxy for phytoplankton biomass (Fig. 3A). For total Chl-a analysis, seawater samples (0.25–1 L) were gently vacuum filtered onto 0.7 µm (nominal porosity) glass fiber filters, which were then stored at −20 °C until analysis. Chl-a concentrations were determined using the acetone extraction and acidification method24,25. Acidification of samples decreased the likelihood of overestimation of Chl-a concentrations due to the presence of chlorophyll degradation products26. Filters were submerged in 10 mL of 90% acetone, sonicated for 10 minutes in an ice bath, and left to extract at −20 °C for 22 h. Following the extraction period, samples were allowed to equilibrate to room temperature (~2 h). Fluorescence of the acetone solution containing the extracted Chl-a was measured before and after acidification with 1.2 N hydrochloric acid using a Turner 10-AU fluorometer. The final concentrations of total Chl-a were calculated from measurements made before (Fo) and after (Fa) acidification using Eq. (1)25. The coefficient (τ) of Eq. (1), adapted from Strickland and Parsons25, was derived from a calibration of the Turner 10-AU fluorometer with known pure chlorophyll standards (Table 2).$${rm{C}}{rm{h}}{rm{l}} mbox{-} {rm{a}}(mu g,{L}^{-1})=frac{tau }{tau -1}ast ({rm{F}}{rm{o}}-{rm{F}}{rm{a}})ast 0.814ast left(frac{{rm{V}}{rm{o}}{rm{l}}.{rm{A}}{rm{c}}{rm{e}}{rm{t}}{rm{o}}{rm{n}}{rm{e}},{rm{e}}{rm{x}}{rm{t}}{rm{r}}{rm{a}}{rm{c}}{rm{t}}{rm{e}}{rm{d}}}{{rm{V}}{rm{o}}{rm{l}}.{rm{S}}{rm{e}}{rm{a}}{rm{w}}{rm{a}}{rm{t}}{rm{e}}{rm{r}},{rm{f}}{rm{i}}{rm{l}}{rm{t}}{rm{e}}{rm{r}}{rm{e}}{rm{d}}}right)$$
(1)
Fig. 3Biological particulate concentrations in the euphotic zone of Saanich Inlet from March 11, 2010 – November 15, 2017. Left panels show time-series profiles (filled contour/scatter plots) of (A) total chlorophyll-a (Total Chl-a), (B) particulate carbon (PC), (C) particulate nitrogen (PN), and (D) particulate biogenic silica (bSiO2). The 2012–2013 data are not interpolated due to single-depth sampling. In A, the bar plot in the top panel shows percent contribution of different size fractions to total Chl-a. In (B–D), black bars in top panels show depth-integrated concentrations. Grey shaded regions in bar plots indicate phytoplankton growing seasons considered for this study (March 1st – October 30th). Right panels show monthly-averaged depth profiles for the entire 8-year period, illustrating euphotic zone seasonality for each particulate. The color scale bars on the far right apply to both the time-series vertical profiles and the 8-year seasonal plots. Sampling depths are indicated by round symbols. The year labels are positioned under the tick marks corresponding to January.Full size imageSize fractionated chlorophyll-aTo determine the percent contributions of “pico” (0.7–2 µm), “small nano” (2–5 µm), “large nano” (5–20 µm) and “micro” ( >20 µm) phytoplankton to total Chl-a, seawater samples (0.25–1 L) separate from those used for total Chl-a) were consecutively filtered through 20, 5 and 2 µm polycarbonate membrane filters and 0.7 µm (nominal porosity) glass fiber filters. Between 2013 and 2017, the “pico” and “small nano” size classes were collected as one fraction (0.7–5 µm). Analysis of Chl-a concentrations for each size fraction followed the same procedure outlined for total Chl-a.Particulate carbon and nitrogenParticulate C and N measurements were obtained from seawater samples incubated for carbon (ρC) and nitrate uptake (ρNO3) rates (see section on “Uptake rates of carbon and nitrate” for methodology) (Fig. 3B,C). PC and PN measurements presented in this dataset were taken at the end of ρC and ρNO3 incubations; however, original (‘ambient’) values can be back calculated by subtracting the amount of C and N taken up during the incubation period from the final PC and PN values. The differences between after-incubation PC and PN data and back-calculated ambient values were not significantly different than the measurement error.Particulate biogenic silicaParticulate biogenic silica was used as a proxy for siliceous phytoplankton biomass (Fig. 3D). Seawater samples (0.5–1 L) from each depth were gently vacuum filtered through 0.6 µm polycarbonate membrane filters. Filters were folded and placed in polypropylene centrifuge tubes, dried for 48 h at 60 °C, and then stored in a desiccator at room temperature until analysis. Filters were digested with 4 mL of 0.2 M NaOH for 30–45 min in a water bath at 95 °C27. After digestion, samples were neutralized with 0.1 N HCl and cooled rapidly in an ice bath. Samples were centrifuged to separate out the undissolved lithogenic silica, and colorimetric analysis was performed on the supernatant. The transmittance of the samples, standards, and reverse-order reagent blanks were read at 820 nm using a Beckman DU 530 ultraviolet-visible (UV/Vis) spectrophotometer27,28.Uptake rates of carbon and nitrateSeawater samples (~0.5–1 L) were gently collected into clear polycarbonate bottles. One additional sample was collected from the 100% light depth, into a dark polycarbonate bottle, which did not allow light penetration. After the addition of the isotopic tracers (see below), bottles were placed into an acrylic incubator with constant seawater flow to maintain surface seawater temperature. Three acrylic tubes wrapped in colored and neutral density photo-film (to obtain 50, 15, and 1% of surface PAR) were used to incubate sample bottles under the same in-situ light conditions from which samples were collected. Samples from the 100% light level were placed inside the same acrylic incubator, but outside of the film-covered tubes. A LI-COR® LI-190 Quantum sensor was installed next to the incubator and continuously recorded incoming PAR for the entire incubation period. During sampling in 2010 and 2011, all experiments were performed using a shipboard incubator. For sampling from 2012 onwards, all experiments were done using an incubator on land (University of Victoria Aquatic Facility), which was connected to a seawater system maintained at local surface seawater temperature (approximately 9–12 °C depending on the time of year).Rates of C (ρC) and NO3 (ρNO3) uptake were determined using a stable isotope tracer-technique29,30 (Fig. 4). A single seawater sample from each light depth received a dual spike, with NaH13CO3 (99% 13C purity, Cambridge Isotope Laboratories) for the determination of ρC and Na15NO3 (98 + % 15N purity, Cambridge Isotopes Laboratories) for the determination of ρNO3. Isotope additions were made at approximately 10% of ambient dissolved inorganic carbon (DIC) and NO3− concentrations.Fig. 4Carbon and nitrate uptake rates in the euphotic zone of Saanich Inlet from March 11, 2010 – November 15, 2017. Left panels show depth-integrated rates (black-bars on top) and time-series profiles (filled contour/scatter plots below) of (A) carbon (ρC) and (B) nitrate (ρNO3) uptake rates. In the time-series profiles, 2012–2013 data are not interpolated due to single-depth sampling. Grey shaded regions in depth-integrated plots indicate phytoplankton growing seasons for this study (March 1st – October 30th). Right panels show monthly-averaged depth profiles for the entire 8-year period, illustrating euphotic zone seasonality for carbon and nitrate uptake. The color scale bars on the far right apply to both the total time-series vertical profiles and the 8-year seasonal plots. Sampling depths are indicated by round symbols. The year labels are positioned under the tick marks corresponding to January.Full size imageSpiked seawater samples were incubated for 24 h, except from 2010 to 2013 when the incubation period was 4 to 6 hr. After incubation, the entire sample was gently vacuum filtered onto a combusted 0.7 µm (nominal porosity) glass fibre filter. Filters were dried for 48 h at 60 °C and kept in a desiccator at room temperature until analysis. Filters were packed into pellets and sent to the Stable Isotope Facility at the University of California (UC) Davis for analysis of 13C and 15N enrichment, and total C and N content by continuous flow isotope ratio mass spectrometry and elemental analysis, respectively. For these measurements, UC Davis uses either an Elementar Vario EL Cube or Micro Cube elemental analyzer (Elementar Analysensysteme GmbH, Hanau, Germany) interfaced to either a PDZ Europa 20–20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, UK) or an Isoprime VisION IRMS (Elementar UK Ltd, Cheadle, UK).Carbon and NO3− uptake rates were calculated using Eq. 3 of Hama et al.29, and Eq. 3 and 6 of Dugdale and Wilkerson30, respectively.For samples incubated for less than 24 h, the daily C or NO3− uptake rates (ρX) were calculated using a PAR extrapolation method shown in Eq. (2):$$rho X(mu mol,{L}^{-1}da{y}^{-1})=left(rho X(mu mol,{L}^{-1}h{r}^{-1})divleft(frac{PAR,during,incubation}{Total,Daily,PAR}right)right)ast 24$$
(2)
Additionally, to account for NO3− uptake occurring under no light, ρNO3− was measured in dark bottles and this rate was added to the ρNO3− of each sample incubated for less than 24 h. The ρNO3− DARK was calculated following Eq. (3):$$rho N{O}_{3}DARK(mu mol,{L}^{-1}da{y}^{-1})=left(rho N{O}_{3}DARK(mu mol,{L}^{-1}h{r}^{-1})divleft(frac{Total,Daily,PAR-PAR,during,incubation}{Total,Daily,PAR}right)right)ast 24$$
(3)
PAR data used in Eqs. (2) and (3) came from the LI-COR® LI-190 Quantum sensor that was mounted beside the incubator. The seawater DIC value for each sample was calculated using a regression equation relating water density to DIC for Saanich Inlet31. Ambient NO3− concentrations were measured as described above. More

An experiment in secondary forests in the Democratic Republic of the Congo finds that calcium, an overlooked soil nutrient, is scarcer than phosphorus, and represents a potentially greater limitation on tropical forest growth.Ecology can reveal distributional patterns and dynamics in nature. One approach used is studying the elemental composition of plants, which has been linked to ecological processes such as growth, diversity or water use efficiency. More recently, elemental composition has been detected as a cofactor in governing the carbon sink capacity of plants, and therefore climate change mitigation1,2,3. This discovery has added an extra layer of urgency to the field, which now aims to better understand and predict global change. The study of nitrogen and/or phosphorus has until now received most of the attention of plant ecologists: nitrogen is the most abundant element in dry leaves after hydrogen and carbon, forms the main structure of proteins and is strongly linked to photosynthesis4. Phosphorus represents around one-tenth of nitrogen abundance in leaves and is key in energy storage and nucleic acids. However, although these represent only two of the many chemical elements that are in flux throughout ecosystems, whether others may have a dominant role in ecosystem dynamics is an open question. Writing in Nature Ecology & Evolution, Bauters et al.5 share some evidence to motivate broadening out from the dominant focus on nitrogen and phosphorus in terrestrial ecology, by revealing a crucial limiting role of calcium in the dynamics of tropical forest succession. More

Gabor transformThe Gabor transform has rarely been used as a feature in a landscape classification OBIA approach but has been used in other OBIA processes such as fingerprint enhancement and human iris detection and for data dimensionality reduction24,29,30,31,32,33,34,35. Gabor filters are a bandpass filter applied to an image to identify texture. The different Gabor bandpass filters mathematically model the visual cortical cells of mammalian brains and thus is expected to improve segmentation and classification accuracy when compared to a human delineated and classified image26,27.Samiappan et al.36 compared Gabor filters to other texture features (grey-level co-occurrence matrix, segmentation-based fractal texture analysis, and wavelet texture analysis) within the GEOBIA process, of a wetland, using sub-meter resolution multispectral imagery. These Gabor filters performed comparably, in overall classification accuracy and Kappa coefficients, with other texture features. However, they were still outperformed by all other texture features. This study did not use any other data for analysis for determining the performance of Gabor filters when paired with data sources such as spectral, NDVI, or LiDAR36,37. Wang et al.38 paired a Gabor transformation with a fast Fourier transformation for edge detection on an urban landscape image that contained uniform textures with promising results. Su30 used the textural attributes derived from Gabor filters for classification but had similar results to Samiappan et al.36 where they found that Gabor features were one of the least useful/influential that contributed to the classification of a mostly agricultural landscape.Gabor filters are a Fourier influenced wavelet transformation, or bandpass filter, that identifies texture as intervals in a 2-D Gaussian modulated sinusoidal wave. This modulation differentiates the Gabor transform from the Fourier transform23,26. These Gabor transformed wavelets are parameterized by the angle at which they alter the image and the frequency of the wavelet. Rather than smoothing an image at the cost of losing detail through Fourier transforms or median filters, Gabor transformed images identify the repeated pattern of localized pixels and gives them similar values if they are a part of the same repeated sequence. Gabor features can closely emulate the visual cortex of mammalian brains that utilize texture to identify objects26,27. This is based on the evaluation of neurons associated with the cortical vertex that respond to different images or light profiles39. Marcelja27 identified that cortical cells responded to signals that are localized frequencies of light like what is represented by the Gabor transformations. Within the frequency domain, the Gabor transform can be defined by Eq. (1):$$Gleft(u, v;f, theta right)= {e}^{-frac{{pi }^{2}}{{f}^{2}} ({gamma }^{2}({u}^{{prime}}-f{)}^{2}+{n}^{2}{v}^{{{prime}}2})}$$
(1)

where (f) is the user-determined frequency (or wavelength); (theta) is the user-determined orientation at which the wavelet is applied to the image; (gamma) and (n) are the standard deviations of the Gaussian function in either direction23,38. These parameters define the shape of the band pass filter and determines its effect on one-dimensional signals. Daugman26, created a 2-D application of this filter in Eq. (2);$$gleft(u,vright)= {e}^{-{pi }^{2}/{f}^{2}[{gamma }^{2}{left({u}^{{prime}}-fright)}^{2}+{n}^{2}{{v}^{{prime}}}^{2}]}$$
(2)

where u’ = ucos − vsin θ θ and v’ = usin − vcos θ.In order to implement Gabor filters on multi-band spectral images, we used Matlab’s Gabor feature on the University of Iowa’s Neon high performance computer (HPC)40 which has up to 512 GB of RAM, which was necessary for processing these images. The first implementation of Gabor filters was performed on a 1610 × 687 single band pixel array (a small subset of the study area), a filter bank of 4 orientations and 8 wavelengths, on a 32 GB RAM computer, and took approximately 8 h to complete. Filter banks are a set of Gabor filters with different parameters that is applied to the spectral image and are required to identify different textures with different orientations and frequencies. By lowering the number of wavelengths from 8 to 4 on an 8128 × 8128 single band pixel array on the same machine 32 GB RAM, the processing was reduced to an hour. Using the HPC, this was further reduced to approximately 90 s using the same filter bank. Before implementing on the HPC, the original spectral image was divided into manageable subsets with overlap in order to prevent ‘edge-effect.’ These images were converted to greyscale by averaging values across all three bands33. When wavelengths become too long, they no longer attribute the textural information desired from the image and therefore add unnecessary computing time. The wavelengths that were used for the filter bank were selected as increasing powers of two starting from 2.82842712475 ((24/sqrt{2})) up to the pixel length of the hypotenuse of the input image. From this, we used only 2.82842712475, 7.0710678, 17.6776695, and 44.19417382. The directional orientation was selected as 45° intervals, from 0 to 180: 0, 45, 90, 135. These parameters were based on the reasoning outlined within Jain and Farrokhina25. More directional orientations could have been included but four were used for computational efficiency. The radial frequencies were selected so that they could capture the different texture in the landscape represented by consistent changes in pixels values within each landcover class. When frequencies are too wide or fine of a width they no longer represent the textures of the different landcover classes and thus are not included. This selection of filter bank parameters are similar or the same as other studies that look into the use of Gabor features for OBIA25,30,31.From the different combinations of parameters (four directions and four frequencies) in the Gabor Transform filter bank, sixteen magnitude response images were created from the converted greyscale three band average image. To limit high local variance within the output Gabor texture images, a Gaussian filter was applied. The magnitude response values were normalized across the 16 different bands so that a Principal Component Analysis (PCA) could be applied. The first principal component of the PCA, from these Gabor transformed images, was used for this study since it limits the computation time to process 16 separate Gabor features, in addition to the other data sources, while still retaining the most amount of information from the different Gabor response features. The Gabor band that was used for this study can be viewed in Fig. 2.Figure 2Gabor transformation. Gabor transformed image of study area derived from original image using the first principal component of all gabor outputs using the filter bank parameters. Software: ArcMap (10.x).Full size imageSegmentationFor this study, we used the watershed algorithm for the segmentation of GEOBIA, implemented by ENVI version 5.0 Feature Extraction tool, due to its ubiquitous use within GEOBIA, its ability to create a hierarchy of segmented objects, and support within the literature as a reliable algorithm37,41,39,43. The watershed algorithm can either use a gradient image or intensity image for segmentation. Based on the observed results, this study used the intensity method. The intensity method averages the value of pixels across bands. Scale, a user-defined parameter, is selected to identify the threshold that decides if a given intensity value within the gradient image can be a boundary. This allows the user to decide the size of the objects created. A secondary, user-defined, parameter defines how similar, adjacent, objects need to be before they are combined or merged. The user arbitrarily selects the parameter value based on how it reduces both under and over segmentation. The parameters selected for this study were visually chosen based on a compromise between over and under segmentation relative to the hand demarcated objects.The merging of two separate objects was based on the full lambda schedule where the user selects a merging threshold ({t}_{i, j}) which is defined by Eq. (3):$${t}_{i, j}= frac{frac{left|{O}_{i}right|cdot left|{O}_{j}right|}{left|{O}_{i}right|+ left|{O}_{j}right|}cdot {Vert {u}_{i}-{u}_{j}Vert }^{2}}{mathrm{length}(mathrm{vartheta }left({O}_{i},{O}_{j}right))}$$
(3)

where ({O}_{i}) is the object of the image, (left|{O}_{i}right|) is the area of (i), ({u}_{i}) is the average of object (i), ({u}_{j}) is the average of object (j), (Vert {u}_{i}-{u}_{j}Vert) is the Euclidean distance between the average values of the pixel values in regions (i) and (j), and (mathrm{length}left(mathrm{vartheta }left({O}_{i},{O}_{j}right)right)) is the length of the shared boundary of ({O}_{i}) and ({O}_{j}).To compare the segmentation of a riparian landscape, with and without Gabor features, we conducted segmentation on two separate sets of data. One dataset was a normalized stacked layer of NDVI and CHM (see Fig. 3) with the original multispectral image used as ancillary data; the other dataset differed only by the inclusion of the Gabor feature. For both instances, the bands were converted to an intensity image by averaging across bands rather than being converted into a gradient image for segmentation. The dataset that included the Gabor features had a scale parameter set at 30 with merge settings at 95 and 95.7 for the sub and super-objects, respectively. The dataset that did not include the Gabor features had a scale parameter of 10 with merge settings at 95.6 and 98.5 for the sub and super-objects, respectively. This resulted in the creation of 87,198 and 62,905 segments for the sub and super objects, respectively, that were created when the Gabor feature was included. 191,050 and 51,664 segments were created for the sub and super objects when the Gabor features, respectively, were not included within the segmentation process. As you will see in the next section, these segments also represent the number of training data that will be included within the supervised classification.Figure 3CHM and NDVI. LiDAR derived canopy height model (top) and normalized difference vegetation index derived from original spectral image. Software: ArcMap (10.x).Full size imageTo create a hierarchy of land cover classes, two sets of segmentation parameters needed to be selected for each dataset. One set of parameters would be used for the sub-objects within the hierarchy and the other set would be used to create super-objects. All parameters used the intensity and full lambda schedule algorithms for the watershed method. The only setting that changed between the sub and super-objects, for either dataset, was the merge parameter which helped maintain similar boundaries as much as possible. Despite this, boundaries could moderately change due to the Euclidean distance, between the pixel values of (i) and (j), changing from the merging of objects; causing ({t}_{i, j}) to cross the threshold which results in a new boundary being drawn. A representation of these results can be viewed and visually compared to the hand demarcated objects in Fig. 4.Figure 4Automated and manual segmented comparison. Juxtaposition of hand delineated, sub-objects, and super-objects for segments generated using the Gabor features. Software: ArcMap (10.x).Full size imageTraining dataThe training data, used for this study, is the transfer of class attributes from hand demarcated and classified segments to automatically segmented objects based on the majority overlap of the hand demarcated segments. Experts identified them using two different classification schemes referenced from the General Wetland Vegetation Classification System44. The 7-class scheme within this system identified objects of either being forest, marsh, agriculture, developed, open water, grass/forbs, or sand/mud. The 13-class scheme identified objects of either being agriculture, developed, grass/forbs, open water, road/levee, sand/mud, scrub-shrub, shallow marsh, submerged aquatic vegetation, upland forest, wet forest, wet meadow, and wet shrub. Not every class from the 7-class scheme will have a sub-class (i.e. developed, open water) but some do for example wet and upland forest are sub-objects of the forest class and wet meadow and shallow marsh are sub-objects of marsh. Figure 5 visually illustrates both classification schemes across the study area.Figure 5Hand delineated objects of both scales. Software: ArcMap (10.x).Full size imageENVI’s feature extraction tool calculates several landscape, spectral, and textural metrics. These attributes were used for each random forest classifier. The Gabor and Hierarchical features will be included selectively to be able to compare their contributions to the (out-of-bag) OOB classification errors. When Gabor features are included within the classification, they are computed the same way as the other image bands.Random forestThe random forest classifier was implemented in R using the random forest module45. The number of trees, that were randomly generated, was large enough (n = 250) to where the Strong law of large numbers would take effect as indicated by the decrease in the change of accuracy. The default number of variables randomly sampled as candidates at each split variable (mtry parameter) was the total number of variables divided by 3 for each dataset. R also generates two separate variable indices: mean decrease in accuracy and mean decrease Gini. Mean decrease in accuracy refers to the accuracy change in the random forest when a single variable is left out. This is a practical metric to determine the usefulness of a variable. The Gini index measures the purity change within a dataset when it is split based upon a given variable within a decision tree.The random forest classification accuracy will be based on the OOB error. The random forest algorithm trains numerous decision trees on random subsets of the training set leaving out a number of training samples when training each decision tree. The samples that are left out of each decision tree are then classified by the decision tree that they were not included within during the training step. The OOB error is the average error of each predicted bootstrapped sample across the ensemble of decision trees within the random forest algorithm.Figure 6 illustrates how the Gabor and hierarchal features were included within the classification of the super and sub-objects.Figure 6Classification procedure. Schematic flow chart illustrating how the Gabor and hierarchal features were included within the classification of the super and sub-objects. OOB classification error included in parenthesis.Full size imageHierarchical schemeTo attribute the hierarchical structure to the sub-objects, we first classified the larger segments that were created with and without the Gabor features using the broader 7-class scheme. These classified super objects were then converted to raster to calculate the majority overlap with the smaller sub-objects. This gave the sub-objects an attribute, the broader 7-class scheme, that could be used to contribute to the classification of the sub-objects with the finer 13-class scheme. This builds the hierarchical relationship between the two class schemes into the supervised classification of the sub-objects. Figure 6 illustrates how the hierarchal structure was included within two of the four sub-object’s list of features used within classification. This methodological approach aligns with O’Neill et al.21 landscape ecology principle that a super-object’s class could be a useful property in defining or predicting a sub-object. This is also different than the more common rule-based approach of iteratively classifying the landscape into smaller and smaller sub-classes22.Segmentation assessmentMost studies rely upon the accuracy assessment of their classifiers to provide support for their analysis results. However, this does not provide evidence whether a new data fusion technique improves the ability to delineate objects of interest within an image. To assess the performance of our segmented polygons, this study evaluated the segments created with and without the Gabor feature using a method highlighted in Xiao et al.37.Our segmentation results were evaluated using an empirical discrepancy measure, used frequently in image segmentation evaluation37,46,47. Discrepancy measures utilize ground truth images that represent the “correct” delineated/classified image to compare the semi-automated image results. In our study, the objects that were delineated and classified by experts from the U.S. Fish and Wildlife Service, were used as training data for our random forest classifier and as ground truth for the discrepancy measure. The discrepancy measure used the percentage of right segmented pixels (PR) in the whole image. To calculate PR, we converted the classified segmented and ground truth polygons to raster and measured the ratio of incorrect pixels to total amount of pixels which was converted to a percentage.Additionally, landscape metrics were calculated using FRAGSTATS48, an open source program commonly used for calculating landscape metrics. FRAGSTATS computed these metrics from thematic raster maps that represent the land cover types of interest. These thematic classes, used for analysis, were the classified objects at both the super and sub-object level. Since we are not attempting to compare the segmentation results for any specific class or area, we calculated metrics on a landscape level. Landscape metrics will represent the segmentation patterns for the entire study area.FRAGSTATS can calculate various metrics representing different aspects of the landscape. The metrics for analysis attempts to understand object geometry. The metrics calculated, for these analyses, were the average and standard deviation for the area (AREA), the fractal dimension index (FRAC), and the perimeter area ratio (PARA). The number of patches (NP) was also included in each result. To take a more landscape centric approach, the area weighted mean was chosen over a simple average. More

Virtual hydrostatic model parametersVarious morphological characteristics were held constant in order to isolate and manipulate the variable of conch shape. A CT-scanned Nautilus pompilius conch was essentially morphed into ammonoid-like conch shapes, populating the Westermann morphospace22 while holding constant septal morphology, septal spacing, and shell/septal thicknesses (Fig. 9). Furthermore, body chamber proportions were determined by iteratively computing soft body volumes that yield Nautilus-like chamber liquid (~ 12% of the phragmocone volume retained)67,68. Septal spacing was measured as the angle from the ventral attachment of the current and previous septa, and the spiraling axis of the conch. Because septal spacing differs in early ontogeny (Fig. S11), only measurements from the 7th to 33rd (terminal) septum were considered. The average angle of 23.46° ± 3.32° (standard deviation) was rounded to 23° and held constant throughout the ontogeny of the hydrostatic models.Figure 9Hydrostatic models of theoretical planispiral cephalopods. These models were constructed by morphing a Nautilus pompilius conch into ammonoid shapes (see “Methods”): (a) oxycone, (b) serpenticone, (c) sphaerocone, and (d) morphospace center. The centers of buoyancy and mass are denoted by the tips of the blue (upper) and red (lower) cones. Prime symbols (′) refer to transparent, transverse views of each respective conch shape. (e) Westermann morphospace22 showing relative positions of these conch shapes. All models were rendered in MeshLab76.Full size imageShell and septal thicknesses were measured with digital calipers from a physical specimen of Nautilus pompilius (Table S13). These measurements were recorded as a ratio of inner whorl height (measured from the ventral point on the current whorl to the ventral point on the previous whorl). These ratios were used in the theoretical models to define shell and septum thicknesses (3.1% of inner whorl height for shell thickness and 2.1% of inner whorl height for septal thickness; Table S13).Hydrostatic model constructionThe near-endmember models were constructed from representative ammonoid specimens (Sphenodiscus lobatus and S. lenticularis—oxycone; Dactylioceras commune—serpenticone; Goniatites crenistria—sphaerocone). Lateral and transverse views were measured from figured specimens for the oxycone (Fig. 5 of Kennedy et al.69), serpenticone (Fig. 2 of Kutygin and Knyazev70), and sphaerocone (Figs. 17 and 20 of Korn and Ebbighausen71). These models were constructed with array algorithms similar to earlier hydrostatic models9,35,72, which were used in a piecewise manner to account for allometric changes in coiling throughout ontogeny (Table S14). These arrays replicated the adult whorl section backwards and translated, rotated, and scaled each successive one. These whorl sections were bridged together to create a single tessellated surface representing the outer interface of the shell. Shell thickness was defined by shrinking the original whorl section so that the thickness between the two was equal to 3.1% of the inner whorl height (Table S13), then using the same array to build the internal interface of the shell. The morphospace center was constructed from previously used conch measurements18 and averaging the whorl section shape in blender (Fig. S12). The corresponding Westermann morphospace parameters (Fig. S13) for each morphology are reported in Table S15.Virtual models of the septa were derived from the CT-scan of Nautilus pompilius (Fig. S14). A single septum was isolated from the adult portion of the phragmocone then smoothed to delete the siphuncular foramen. This septum was placed within the whorl section of each theoretical model and stretched in the lateral directions until it approximately fit. The “magnetize” tool in Meshmixer (Autodesk Inc.) was used to attach the septal margin to the new whorl section so that the Nautilus suture was transferred to the new whorl section. The septum was then smoothed to reconcile the first order curves with the new location of the septal margin. The respective septum for each theoretical model was then replicated with the same array instructions used to build the shell. Because each replicated object was rotated one degree (Table S14), 22 septa were deleted in between every two so that the septal spacing was equal to 23° (Fig. S11).For each theoretical model, the septa were unified with the model of the shell using Boolean operations in Netfabb (Autodesk Inc.). To perform hydrostatic calculations, virtual models must be created for each material of unique density. The virtual model of the shell constrains the shape of the soft body (within the body chamber) and chamber volumes (within the phragmocone). These internal interfaces were isolated from the model of the shell, then their faces inverted for proper, outward-facing orientations of their normals. A conservative soft body estimate was created, aligning with previously published reconstructions64,65,73. The profile shape of this soft body was scaled and maintained between each model. External interfaces of the shell and soft body were also isolated to create a model of the water displaced by each theoretical cephalopod. Each of these models are necessary for hydrostatic calculations (buoyancy and the distribution of organismal mass).Each hydrostatic model is stored in an online repository (Dataset S1; https://doi.org/10.5281/zenodo.5684906). The hydrostatic centers of each virtual model and their volumes and masses are listed in Tables S16 and S17.Hydrostatic calculationsEach theoretical model was scaled to have equal volume (near one kilogram; 0.982 kg–a result of arbitrarily scaling the sphaerocone model to 15 cm in conch diameter). An object is neutrally buoyant when the sum of organismal mass is equal to the mass of water displaced (the principle of Archimedes). The percentage of chamber liquid can be computed to satisfy this condition.$${Phi } = frac{{left( {frac{{{text{V}}_{{{text{wd}}}} {uprho }_{{{text{wd}}}} – {text{V}}_{{{text{sb}}}} {uprho }_{{{text{sb}}}} – {text{V}}_{{{text{sh}}}} {uprho }_{{{text{sh}}}} }}{{{text{V}}_{{{text{ct}}}} }}} right) – left( {{uprho }_{{{text{cl}}}} } right)}}{{left( {{uprho }_{{{text{cg}}}} – {uprho }_{{{text{cl}}}} } right)}}$$
(1)
where Vwd and ρwd are the volume and density of the water displaced, Vsb and ρsb are the volume and density of the soft body, Vsh and ρsh are the volume and density of the shell, ρcl is the density of cameral liquid, ρcg is the density of cameral gas, and Vct is the total volume of all chambers. A soft body density of 1.049 g/cm3 is used based on bulk density calculations of Nautilus-like tissues74, a seawater-filled mantle cavity, and thin calcitic mouthparts21. A shell density of 2.54 g/cm374, cameral liquid density of 1.025 g/cm375, and cameral gas density of 0.001 g/cm3 are adopted from recent hydrostatic studies.Other hydrostatic properties depend on the relative positions of the centers of buoyancy and mass. The center of buoyancy is equal to the center of volume of water displaced. This center and the centers of each virtual model of unique density were computed in the program MeshLab76. The individual centers for each organismal model (soft body, shell, cameral liquid and cameral gas) were used to compute the total center of mass, with an average weighted by material density:$$M = frac{{sum left( {L*m_{o} } right)}}{{sum m_{o} }}$$
(2)
where M is the total center of mass in a principal direction, L is the center of mass of a single object measured with respect to an arbitrary datum in each principal direction, and (m_{o}) is the mass of each object with unique density. Equation 2 was used in the x, y, and z directions to compute the 3D coordinate position of the center of mass. The centers of mass for the chamber contents (liquid and gas) were set equal to the center of volume of all chambers, a minor assumption considering the capillary retention of liquid around the septal margins in the living animals62.The hydrostatic stability index (St) is computed from the relative location of the centers of buoyancy (B) and mass (M), normalized by the cube root of volume (V) for a dimensionless metric that is independent of scale:$$S_{t} = frac{{ sqrt {left( {B_{x} – M_{x} } right)^{2} + left( {B_{y} – M_{y} } right)^{2} + left( {B_{z} – M_{z} } right)^{2} } }}{{sqrt[3]{V}}}$$
(3)
where the subscripts correspond to the x, y, and z components of each hydrostatic center.Apertural orientations were measured in blender after orienting each model so that the center of buoyancy was vertically aligned above the center of mass. Apertural angles of 0° correspond to a horizontally facing soft body, while angles of + 90° and − 90° correspond to upward- and downward-facing orientations, respectively.Thrust angles were measured from the hyponome location (ventral edge of the aperture) to the midpoint of the hydrostatic centers, with respect to the horizontal. Thrust angles of 0° infer idealized horizontal backward transmission of energy into movement, while thrust angles of + 90° and − 90° infer more efficient transmission of energy into downward and upward vertical movement, respectively.Biomimetic robot constructionTo isolate the variable of shell shape on swimming capabilities, only the external shape, and static orientation of each virtual hydrostatic model were used to build physical, 3D printed robots. That is, each model has artificially high hydrostatic stability (Tables S3) to nullify the effect of the thrust angle (the angle at which thrust energy passes through the hydrostatic centers and most efficiently transmits energy into movement; Table S4). Less stable morphotypes (e.g., serpenticones and sphaerocones) are more sensitive to the constraints imposed by this hydrostatic property.Space constraints inside each model were determined by first constructing a propulsion system and electronic components that operate the motor. The models use impeller-based water pumps (Figs. 1d and 10a) driven by a brushed DC motor. This system creates a partial vacuum by centrifugal acceleration, drawing water from a “mantle cavity” and expelling it out of a “hyponome”. This system was optimized by iteratively designing models in Blender77, then testing 3D-printed, stand-alone water pumps. After three iterations, a four-blade impeller and gently tapering hyponome (inner diameter at distal end = 6.7 mm) were chosen. The electronic components used to drive the motor consist of an Arduino Pro Micro microcontroller, a motor driver, and two batteries (Fig. 10). A 3.7 V battery operates the microcontroller, and a larger 7.4 V battery supplies power to the motor. Communication is achieved via infrared, allowing specification of the jet pulse duration, number of pulses, and the power level of the motor (using pulse-width modulation; PWM). Each of these electronic components fold into a compact cartridge capable of being plugged into 3D-printed models of each investigated shell shape (Figs. 2 and 10). Each model was designed with brackets to hold the electronics cartridge in place. The sphaerocone had the most severe space constraints, with low conch diameter to volume ratio. After determining the space required for the electronics (Fig. 10) this model was scaled to 15 cm, and all other models were scaled to have similar volumes (with subtle volume differences due to minor differences in soft body shape compared to the hydrostatic models).Figure 10Biomimetic cephalopod robot components. (a) Ventral view of the sphaerocone biomimetic robot (before covering the pump and mantle cavities) with assembled electronics cartridge to the right. (b) View of electronic components that fit into the cartridge. (c) Electronics cartridge placed in robot. These two halves are fit together with wax to create a water-tight seal. Each model component is denoted by letters in circles: A = Arduino microcontroller, B = microcontroller charger / voltage regulator, C = motor driver, D = infrared sensor, E = indicator LED, F = microcontroller battery (3.7 V), G = motor battery (7.4 V), H = brushed motor, I = impeller and water pump cavity, J = electronics cartridge. The colors of annotations correspond to components depicted in Figs. 1 and 2.Full size imageIn addition to having a propulsion system, biomimetic cephalopod robots must also be capable of neutral buoyancy, while assuming the proper orientation in the water. These robots, and their once-living counterparts, each have differing material densities and associated mass distributions for each component. To reconcile these differences, the total mass and total centers of mass for each model were manipulated by controlling the volume and 3D distribution of the 3D-printed PETG (polyethylene terephthalate glycol) thermoplastic. That is, the shape of this material holds each model component in place while correcting for these differences in hydrostatics. The PETG mass required for neutral buoyancy was found by subtracting the mass of every other model component from the mass of the water displaced by the model (i.e., electronics cartridge, bismuth counterweight, liquid, motor, batteries, electronic components, and self-healing rubber; Table S1). This model configuration also allows buoyancy to be fine-tuned in water, compensating for potential density differences between the virtual water and the actual water in the experimental settings. That is, each virtual model accounts for ~ 9 g of internal liquid, but the actual volume of this liquid can be adjusted in the physical robot with a syringe through a self-healing rubber valve (Table S1; Fig. 1).The 3D position of the total center of mass was manipulated by accounting for the local centers of mass of each material of unique density. Materials like the batteries, motor, and electronic components were each assigned bulk density values because they are made up of composite materials. While this is an approximation, their contributions to the total center of mass are low because they account for small fractions of the total model mass (Tables S1 and S2). These components, like all others, were digitally modeled in Blender77 and their volumes and centers of mass were computed in the program MeshLab76. A dense, bismuth counterweight was also modeled, and positioned to artificially stabilize each model (pulling the z component of the total center of mass downward, while maintaining the horizontal components). The virtual model of this counterweight was used to make a 3D-printed mold, allowing a high heat silicone mold to be casted. The bismuth counterweight was cast from this silicone mold and filed/sanded to the dimensions of its virtual counterpart. Hyponomes were oriented horizontally, to yield movement in this direction. To maintain the same static orientation as the virtual model (same x and y center of mass components), the PETG center of mass was computed with the following equation:$$D_{PETG} = frac{{Mmathop sum nolimits_{i = 1}^{n} m_{i} – mathop sum nolimits_{i = 1}^{n} (D_{i} m_{i} )}}{{left( {m_{PETG} } right)}}$$
(4)
where DPETG is the location of the PETG center of mass from an arbitrary datum in each principal direction. M is the total center of mass in a particular principal direction, mi is the mass of each model component, Di is the local center of mass of each model component in a particular principal direction and mPETG is the mass of the PETG required for a neutrally buoyant condition. See Tables S1 and S2 for a list of model components and measurements.Each model was 3D printed with an Ultimaker S5 3D printer using clear (natural) PETG in separate parts, allowing the internal components to be implanted (i.e., brushed DC motors and bismuth counterweights). Each model part was chemically welded together with 100% dichloromethane, with minor amounts of cyanoacrylate glue used to fill seams (e.g., the water pump lid; Fig. 10a). Each final model consists of the main body (housing the water pump, motor, and counterweight), and a “lid” with brackets that house the electronics cartridge (Figs. 2 and 10). The main body and lid were fused together before each experiment by placing wax (paraffin-beeswax blend) along a tongue and groove seam, heating it with a hairdryer, then vigorously squeezing each part together. Surplus wax extruded from the seam was removed and smoothed, producing a water-tight seal.Thrust calibrationEven though each model was designed to have equal mantle cavity and pump cavity volumes, they produced slightly different thrusts. These differences were likely due to variable degrees of friction between the impellers and the surrounding water pumps. To correct for these differences, the thrust produced by each model was measured with a Vernier Dual-Range Force Sensor (0.01 N resolution). Each robot was attached at the hyponome location, through a series of pulleys, and to the sensor with fishing line (Fig. S1; similar to the methods used for living cephalopods78). Force was recorded for 30-s intervals at a sample rate of 0.05 s. During this time, each model was recorded jetting with a 6-s pulse for 15 trials (Fig. S2A). Each trial had initial noise from setting up the model, then peaked randomly when the fishing line became taught, then stabilized after some period of oscillation. Only the stabilized portion of the thrust profile was used to record thrust at 100% voltage for each model (Fig. S2B). The true zero datum was also subtracted from each of these trials. The lowest thrust from each of the models was used as a baseline (serpenticone and oxycone). Each model was recorded again for 15 trials by lowering the motor voltage in increments of 5% until they yielded similar thrusts (0.3 N) to the original serpenticone and oxycone trials (Fig. S2C). The final power levels were then determined for each model and adjusted with pulse-width modulation (PWM) through the microcontroller: serpenticone (100%), oxycone (100%), sphaerocone (95%), and morphospace center (85%).The peak thrust measured for 1 kg extant Nautilus is around 2 N16. The time-averaged thrust during each pulse is around 23% of this value (0.46 N16). This computed value slightly overpredicts observed maximum velocities for this animal (33 cm/s instead of 25 cm/s), so the appropriate time-averaged thrust is probably slightly lower. The motor in the robots quickly reaches its maximum thrust (~ 0.3 N) once initiated then quickly declines after shutting off (Fig. S2). Therefore, the thrust produced by the robots can be treated as a conservative Nautilus-like jet thrust close to the behavior of escape jetting. One-second pulse and refill intervals are also on par with values reported for extant Nautilus16.Robot buoyancyEach of the models were made near neutrally buoyant by adjusting the allotted ~ 9 g of internal liquid with a syringe through a self-healing rubber valve. The single-pulse experiments were performed in an external pool (ranging ~ 23.5 to 26.5 °C). The three-pulse and maneuverability experiments were performed in an internal pool (the Crimson Lagoon at the University of Utah). This internal pool had slightly higher temperatures (~ 28 °C), yielding lower ambient water densities than the virtual water. These conditions required slightly less internal liquid (~ 2–5 g). These differences in internal liquid masses produced negligibly small shifts in mass distributions because they are very small proportions of total robot masses (Table S1).Perfect neutral buoyancy cannot be practically achieved, but this condition can be closely approached. Each of the biomimetic robots experience subtle upward or downward movements of the course of their 5–15 s long trials due to slightly positive or negative buoyancies. Because these differences in buoyancy influence the vertical component of movement, only the horizontal components are considered for discussion. However, a comparison of velocities computed from full, 3D movement (Eq. 5) and restricted 2D components (Eq. 6) reveals that these differences are minor (Figs. S7 and S8). These comparisons demonstrate that model buoyancy did not substantially influence kinematics other than gross trajectories (Figs. 4 and S9).3D motion trackingAfter adjusting buoyancy, each model was positioned underwater with a grabber tool. This tool was fitted with a bundle of fiber-optic cable (Fig. S4) attached to an infrared remote control. Arduino code (Dataset S2) was uploaded to the microcontroller in the robot allowing jet pulse duration, number of pulses, and power to be adjusted with this remote control. After an infrared pulse is received, the motor activates, and activity is indicated by a green LED that illuminates the model from the inside. This light is used to determine time-zero for each trial of motion tracking.After sending an infrared signal, the movement of each model was recorded with a submersible camera rig fitted with two waterproof cameras (Fig. 3). Each of the four models were monitored during a single, one-second jet for at least 9 trials each. Additionally, the laterally compressed morphotypes (serpenticone and oxycone) were monitored during three, one-second pulses for 10 trials each. The inflated morphotypes (sphaerocone and morphospace center) were not able to be monitored over longer distances because they had the tendency to rotate about the vertical axis, obscuring views of the tracking points. In addition to horizontal movement, turning efficiency (maneuverability about the vertical axis) was monitored by directing the cameras with a top-down view of each model. A 90° elbow attachment for the hyponome was fit to each model to investigate the ease or difficulty of rotation. Each model was designed to spin counter-clockwise when viewed from above so that the influence of the motor’s angular momentum was consistent between models.Footage was recorded with two GoPro Hero 8 Black cameras at 4K resolution and 24 (23.975) frames per second, with linear fields of view. Motion tracking was performed with the software DLTdv879 to record the pixel locations of each tracking point (Figs. 1c and S4). These coordinates were transformed into 3D coordinates in meters using the program easyWand580. The tracking points on each model were used for wand calibration because the distances between these sets of points were fixed. Standard deviations of the reproduced tracking point distances of less than 1 cm were considered suitable.The 3D position datasets allowed velocity, acceleration, rocking, to be computed for each experiment. Additionally angular displacement and angular velocity was of interest for the rotation experiments about the vertical axis. Velocity was computed under two scenarios: (1) using the 3D movement direction between each timestep (Eq. 5), and (2) only considering the horizontal movement direction between each time step (Eq. 6). The latter scenario was preferred to nullify the influences of model buoyancies, which were not perfectly neutral and caused some degree of vertical movement.$$V_{i} = frac{{sqrt {left( {x_{i} – x_{i – 1} } right)^{2} + left( {y_{i} – y_{i – 1} } right)^{2} + left( {z_{i} – z_{i – 1} } right)^{2} } }}{{left( {t_{i} – t_{i – 1} } right)}}$$
(5)
$$V_{i} = frac{{sqrt {left( {x_{i} – x_{i – 1} } right)^{2} + left( {y_{i} – y_{i – 1} } right)^{2} } }}{{left( {t_{i} – t_{i – 1} } right)}}$$
(6)
where V and t are velocity and time, and the subscripts i and i −1 refer to the current and previous time steps, respectively. Coordinate components are denoted by x, y, and z at each timestep. The averaged 3D location of both tracking points was used for each model (i.e., midpoints). Note that Eq. (5) uses the 3D form of the Theorem of Pythagoras, whereas Eq. (6) uses the 2D version. Time zero for each trial was defined as the frame where the robot was illuminated by the internal LED, indicating motor activity. Acceleration was modeled by fitting a linear equation to the datapoints during the one-second pulse interval(s) using the curve fitting toolbox in MATLAB R2020A.The artificially high hydrostatic stability of each model was designed to nullify rocking during movement. This behavior was computed for each model during the one-pulse experiments with the following equation:$$theta_{dv} = cos^{ – 1} left( {frac{{left( {z_{2} – z_{1} } right)}}{{sqrt {left( {x_{2} – x_{1} } right)^{2} + left( {y_{2} – y_{1} } right)^{2} + left( {z_{2} – z_{1} } right)^{2} } }}} right) – theta_{tp}$$
(7)
where (theta_{dv}) is the angle deviated from true vertical and (theta_{tp}) is the angle of the tracking points measured from the vertical in a static setting. The subscripts 1 and 2 of the x, y, and z coordinates refer to the anterior and posterior tracking points, respectively.Maneuverability about the vertical axis was determined by computing the angle between the horizontal components of each tracking point. The net angle from the starting angle for each trial was tabulated. Angular velocity was determined by dividing the change in angle between each frame by the frame duration (1/23.975 fps).Links to example motion tracking footage, and robotic models are deposited in an online repository60,61,63 (Dataset S2; https://doi.org/10.5281/zenodo.6180801). More

Sabkha soil analysisSoil biochemical analysis shows higher concentration of salt in the middle (average of 200 ppt) which decreases to 20 ppt outside of the sabkha (Fig. 1B). Soil conductivity and total dissolved solids (TDS) were positively correlated with salinity (r = 0.8, p  More