Ecology

Kim, S. Y. & Yi, D. Y. Components of human breast milk: From macronutrient to microbiome and microRNA. Clin. Exp. Pediatr. 63(8), 301 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Power, M. L. & Schulkin, J. Maternal regulation of offspring development in mammals is an ancient adaptation tied to lactation. Appl. Transl. Genomics. 2, 55–63 (2013).CAS 
Article 

Google Scholar 
Pannaraj, P. S. et al. Association between breast milk bacterial communities and establishment and development of the infant gut microbiome. JAMA Pediatr. 171(7), 647–654 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Lyons, K. E., Ryan, C. A., Dempsey, E. M., Ross, R. P. & Stanton, C. Breast milk, a source of beneficial microbes and associated benefits for infant health. Nutrients 12(4), 1039 (2020).CAS 
PubMed Central 
Article 

Google Scholar 
Fehr, K. et al. Breastmilk feeding practices are associated with the co-occurrence of bacteria in mothers’ milk and the infant gut: The CHILD cohort study. Cell Host Microbe. 28(2), 285–297 (2020).CAS 
PubMed 
Article 

Google Scholar 
Moossavi, S. & Azad, M. B. Origins of human milk microbiota: New evidence and arising questions. Gut Microbes. 12(1), 1667722. https://doi.org/10.1080/19490976.2019.1667722 (2020).CAS 
Article 
PubMed 

Google Scholar 
Groer, M. W., Morgan, K. H., Louis-Jacques, A. & Miller, E. M. A scoping review of research on the human milk microbiome. J. Hum. Lact. 36(4), 628–643 (2020).PubMed 
Article 

Google Scholar 
Gopalakrishna, K. P. & Hand, T. W. Influence of maternal milk on the neonatal intestinal microbiome. Nutrients 12(3), 823 (2020).CAS 
PubMed Central 
Article 

Google Scholar 
Ayoub Moubareck, C., Lootah, M., Tahlak, M. & Venema, K. Profiles of human milk oligosaccharides and their relations to the milk microbiota of breastfeeding mothers in Dubai. Nutrients 12(6), 1727 (2020).PubMed Central 
Article 
CAS 

Google Scholar 
Henrick, B. M. et al. Bifidobacteria-mediated immune system imprinting early in life. Cell 184, 3884–3898 (2021).CAS 
PubMed 
Article 

Google Scholar 
Walker, W. A. & Iyengar, R. S. Breast milk, microbiota, and intestinal immune homeostasis. Pediatr. Res. 77(1), 220–228 (2015).CAS 
PubMed 
Article 

Google Scholar 
Petrullo, L. et al. The early life microbiota mediates maternal effects on offspring growth in a nonhuman primate. Iscience. 25(3), 103948 (2022).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bowen, W. D., Boness, D. J. & Oftedal, O. T. Mass transfer from mother to pup and subsequent mass loss by the weaned pup in the hooded seal, Cystophora cristata. Can. J. Zool. 65(1), 1–8 (1987).Article 

Google Scholar 
Smith, T. M., Austin, C., Hinde, K., Vogel, E. R. & Arora, M. Cyclical nursing patterns in wild orangutans. Sci. Adv. 3(5), e1601517 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Park, Y. W. & Haenlein, G. F. W. Handbook of Milk of Non-Bovine Mammals (Wiley, 2008).
Google Scholar 
Oftedal, O. T. Use of maternal reserves as a lactation strategy in large mammals. Proc. Nutr. Soc. 59(1), 99–106 (2000).CAS 
PubMed 
Article 

Google Scholar 
Hinde, K. & Milligan, L. A. Primate milk: Proximate mechanisms and ultimate perspectives. Evol. Anthropol. Issues News Rev. 20(1), 9–23 (2011).Article 

Google Scholar 
Osthoff, G., Hugo, A., De Wit, M., Nguyen, T. P. M. & Seier, J. Milk composition of captive vervet monkey (Chlorocebus pygerythrus) and rhesus macaque (Macaca mulatta) with observations on gorilla (Gorilla gorilla gorilla) and white handed gibbon (Hylobates lar). Comp. Biochem. Physiol. Part B Biochem. Mol. Biol. 152(4), 332–338 (2009).CAS 
Article 

Google Scholar 
Power, M. L., Oftedal, O. T. & Tardif, S. D. Does the milk of callitrichid monkeys differ from that of larger anthropoids?. Am. J. Primatol. Off. J. Am. Soc. Primatol. 56(2), 117–127 (2002).
Google Scholar 
Power, M. L. et al. Patterns of milk macronutrients and bioactive molecules across lactation in a western lowland gorilla (Gorilla gorilla) and a Sumatran orangutan (Pongo abelii). Am. J. Primatol. 79(3), e22609 (2017).Article 
CAS 

Google Scholar 
Garcia, M., Power, M. L. & Moyes, K. M. Immunoglobulin A and nutrients in milk from great apes throughout lactation. Am. J. Primatol. 79(3), e22614 (2017).Article 
CAS 

Google Scholar 
Muletz-Wolz, C. R. et al. Diversity and temporal dynamics of primate milk microbiomes. Am. J. Primatol. 81(10–11), e22994 (2019).PubMed 
PubMed Central 

Google Scholar 
Rodríguez, J. M. The origin of human milk bacteria: Is there a bacterial entero-mammary pathway during late pregnancy and lactation?. Adv. Nutr. 5(6), 779–784 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
LaTuga MS, Stuebe A, Seed PC. A review of the source and function of microbiota in breast milk. In Seminars in Reproductive Medicine, Vol 32, 68–73 (Thieme Medical Publishers, 2014).Chen, W. et al. Lactation stage-dependency of the sow milk microbiota. Front. Microbiol. 9, 945 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
McInnis, E. A., Kalanetra, K. M., Mills, D. A. & Maga, E. A. Analysis of raw goat milk microbiota: Impact of stage of lactation and lysozyme on microbial diversity. Food Microbiol. 46, 121–131 (2015).CAS 
PubMed 
Article 

Google Scholar 
Gonzalez, E. et al. Distinct changes occur in the human breast milk microbiome between early and established lactation in breastfeeding Guatemalan mothers. Front. Microbiol. 12, 194 (2021).Article 

Google Scholar 
Ge, Y. et al. The maternal milk microbiome in mammals of different types and its potential role in the neonatal gut microbiota composition. Animals 11(12), 3349 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Kordy, K. et al. Contributions to human breast milk microbiome and enteromammary transfer of Bifidobacterium breve. PLoS ONE 15(1), e0219633 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Jost, T., Lacroix, C., Braegger, C. & Chassard, C. Impact of human milk bacteria and oligosaccharides on neonatal gut microbiota establishment and gut health. Nutr. Rev. 73(7), 426–437 (2015).PubMed 
Article 

Google Scholar 
Fernández, L. et al. The human milk microbiota: Origin and potential roles in health and disease. Pharmacol. Res. 69(1), 1–10 (2013).PubMed 
Article 
CAS 

Google Scholar 
Cabrera-Rubio, R. et al. The human milk microbiome changes over lactation and is shaped by maternal weight and mode of delivery. Am. J. Clin. Nutr. 96(3), 544–551 (2012).CAS 
PubMed 
Article 

Google Scholar 
Gomez-Gallego, C., Garcia-Mantrana, I., Salminen, S. & Collado, M. C. The human milk microbiome and factors influencing its composition and activity. In Seminars in Fetal and Neonatal Medicine. Vol 21, 400–405 (Elsevier, 2016).Khodayar-Pardo, P., Mira-Pascual, L., Collado, M. C. & Martínez-Costa, C. Impact of lactation stage, gestational age and mode of delivery on breast milk microbiota. J. Perinatol. 34(8), 599–605 (2014).CAS 
PubMed 
Article 

Google Scholar 
Wan, Y. et al. Human milk microbiota development during lactation and its relation to maternal geographic location and gestational hypertensive status. Gut Microbes. 11(5), 1438–1449 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hunt, K. M. et al. Characterization of the diversity and temporal stability of bacterial communities in human milk. PLoS ONE 6(6), e21313 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Petrullo, L., Jorgensen, M. J., Snyder-Mackler, N. & Lu, A. Composition and stability of the vervet monkey milk microbiome. Am. J. Primatol. 81(10–11), e22982 (2019).PubMed 

Google Scholar 
Mittermeier, R. A. et al. Primates in peril: The world’s 25 most endangered primates 2008–2010. Primate Conserv. 24(1), 1–57 (2009).Article 

Google Scholar 
Williams, J. E. et al. Human milk microbial community structure is relatively stable and related to variations in macronutrient and micronutrient intakes in healthy lactating women. J. Nutr. 147(9), 1739–1748 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
Kumar, H. et al. Distinct patterns in human milk microbiota and fatty acid profiles across specific geographic locations. Front. Microbiol. 7, 1619 (2016).PubMed 
PubMed Central 

Google Scholar 
Keady, M. et al. Clinical health issues, reproductive hormones, and metabolic hormones associated with gut microbiome structure in African and Asian elephants. Anim. Microbiome. 3, 1–19 (2021).Article 
CAS 

Google Scholar 
RStudio Team. RStudio: Integrated Development for R. http://www.rstudio.com/ (2020).Bolyen, E. et al. QIIME 2: Reproducible, Interactive, Scalable, and Extensible Microbiome Data Science. PeerJ Preprints (2018).Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods. 13(7), 581 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Cole, J. R. et al. Ribosomal Database Project: Data and tools for high throughput rRNA analysis. Nucleic Acids Res. 42(D1), D633–D642 (2014).CAS 
PubMed 
Article 

Google Scholar 
Davis, N. M., Proctor, D. M., Holmes, S. P., Relman, D. A. & Callahan, B. J. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome. 6(1), 1–14 (2018).Article 

Google Scholar 
McMurdie, P. J. & Holmes, S. phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8(4), e61217 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Beule, L. & Karlovsky, P. Improved normalization of species count data in ecology by scaling with ranked subsampling (SRS): Application to microbial communities. PeerJ 8, e9593 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V. & Egozcue, J. J. Microbiome datasets are compositional: And this is not optional. Front. Microbiol. 8, 2224 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Oksanen, J. et al. vegan: Community Ecology Package. https://cran.r-project.org/package=vegan (2020).Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57, 289–300 (1995).MathSciNet 
MATH 

Google Scholar 
Kumbhare, S. V., Patangia, D. V., Patil, R. H., Shouche, Y. S. & Patil, N. P. Factors influencing the gut microbiome in children: From infancy to childhood. J. Biosci. 44(2), 1–19 (2019).Article 

Google Scholar 
Amato, K. R. et al. Phylogenetic and ecological factors impact the gut microbiota of two Neotropical primate species. Oecologia 180(3), 717–733 (2016).ADS 
PubMed 
Article 

Google Scholar 
Mulligan, M. E. et al. Methicillin-resistant Staphylococcus aureus: A consensus review of the microbiology, pathogenesis, and epidemiology with implications for prevention and management. Am. J. Med. 94(3), 313–328 (1993).CAS 
PubMed 
Article 

Google Scholar 
Ruegg, P. L. A 100-Year Review: Mastitis detection, management, and prevention. J. Dairy Sci. 100(12), 10381–10397 (2017).CAS 
PubMed 
Article 

Google Scholar 
Clarridge, J. E. III. Impact of 16S rRNA gene sequence analysis for identification of bacteria on clinical microbiology and infectious diseases. Clin. Microbiol. Rev. 17(4), 840–862 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Martín, V., Mediano, P., Del Campo, R., Rodríguez, J. M. & Marín, M. Streptococcal diversity of human milk and comparison of different methods for the taxonomic identification of streptococci. J. Hum. Lact. 32(4), NP84–NP94 (2016).PubMed 
Article 

Google Scholar 
Ghebremedhin, B., Layer, F., Konig, W. & Konig, B. Genetic classification and distinguishing of Staphylococcus species based on different partial gap, 16S rRNA, hsp60, rpoB, sodA, and tuf gene sequences. J. Clin. Microbiol. 46(3), 1019–1025 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Chen, Q. et al. Quantification of human oral and fecal Streptococcus parasanguinis by use of quantitative real-time PCR targeting the groEL gene. Front. Microbiol. 10, 2910 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Plows, J. F. et al. Longitudinal changes in human milk oligosaccharides (HMOs) over the course of 24 months of lactation. J. Nutr. 151(4), 876–882 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Boehm, G. & Stahl, B. Oligosaccharides from milk. J. Nutr. 137(3), 847S-849S (2007).CAS 
PubMed 
Article 

Google Scholar 
van Leeuwen, S. S. et al. Goat milk oligosaccharides: Their diversity, quantity, and functional properties in comparison to human milk oligosaccharides. J. Agric. Food Chem. 68(47), 13469–13485 (2020).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Tao, N. et al. Evolutionary glycomics: Characterization of milk oligosaccharides in primates. J. Proteome Res. 10(4), 1548–1557 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Yu, Z.-T., Chen, C. & Newburg, D. S. Utilization of major fucosylated and sialylated human milk oligosaccharides by isolated human gut microbes. Glycobiology 23(11), 1281–1292 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bolotin, A. et al. Complete sequence and comparative genome analysis of the dairy bacterium Streptococcus thermophilus. Nat. Biotechnol. 22(12), 1554–1558 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Schwab, C. & Gänzle, M. Lactic acid bacteria fermentation of human milk oligosaccharide components, human milk oligosaccharides and galactooligosaccharides. FEMS Microbiol. Lett. 315(2), 141–148 (2011).CAS 
PubMed 
Article 

Google Scholar 
Marcobal, A. et al. Consumption of human milk oligosaccharides by gut-related microbes. J. Agric. Food Chem. 58(9), 5334–5340 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Uriot, O. et al. Streptococcus thermophilus: From yogurt starter to a new promising probiotic candidate?. J. Funct. Foods. 37, 74–89 (2017).CAS 
Article 

Google Scholar 
Duar, R. M., Henrick, B. M., Casaburi, G. & Frese, S. A. Integrating the ecosystem services framework to define dysbiosis of the breastfed infant gut: The role of B. infantis and human milk oligosaccharides. Front. Nutr. 7, 33 (2020).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Singh, R. P., Niharika, J., Kondepudi, K. K., Bishnoi, M. & Tingirikari, J. M. R. Recent understanding of human milk oligosaccharides in establishing infant gut microbiome and roles in immune system. Food Res. Int. 151, 110884. https://doi.org/10.1016/j.foodres.2021.110884 (2022).CAS 
Article 
PubMed 

Google Scholar 
Ximenez, C. & Torres, J. Development of microbiota in infants and its role in maturation of gut mucosa and immune system. Arch. Med. Res. 48(8), 666–680. https://doi.org/10.1016/j.arcmed.2017.11.007 (2017).Article 
PubMed 

Google Scholar 
Meehan, C. L. et al. Social networks, cooperative breeding, and the human milk microbiome. Am. J. Hum. Biol. 30(4), e23131 (2018).PubMed 
Article 

Google Scholar 
Bornbusch, S. L. et al. Stable and transient structural variation in lemur vaginal, labial and axillary microbiomes: Patterns by species, body site, ovarian hormones and forest access. FEMS Microbiol. Ecol. 96(6), fiaa090 (2020).CAS 
PubMed 
Article 

Google Scholar 
Bornbusch, S. L. & Drea, C. M. Antibiotic resistance genes in lemur gut and soil microbiota along a gradient of anthropogenic disturbance. Front. Ecol. Evol. https://doi.org/10.3389/fevo.2021.704070 (2021).Article 

Google Scholar 
Grieneisen, L. E. et al. Genes, geology and germs: Gut microbiota across a primate hybrid zone are explained by site soil properties, not host species. Proc. R. Soc. B. 2019(286), 20190431 (1901).
Google Scholar 
Ellison, S. et al. The influence of habitat and phylogeny on the skin microbiome of amphibians in Guatemala and Mexico. Microb. Ecol. 78(1), 257–267 (2019).PubMed 
Article 

Google Scholar 
Phillips, C. D. et al. Microbiome analysis among bats describes influences of host phylogeny, life history, physiology and geography. Mol. Ecol. 21(11), 2617–2627 (2012).PubMed 
Article 

Google Scholar  More

Dynamic in-flight flock organizationIt is commonly assumed that during flocking, flock members follow three basic interaction rules: Attraction, Repulsion and Alignment, to coordinate spatial positions between each other18. To study the spatial organization of our zebra finch flock during flight, the spatial positions of all birds in the flight section were tracked in every fifth frame (sample rate: 24 Hz (that is, frames per second)) of the synchronized footage recorded by two high-speed digital video cameras (Camera 1: centred upwind view, Fig. 1a,b; Camera 2: upturned vertical view, Fig. 1a,c) for the entire duration (51.7, 58.3, 69.2 and 127 s) of four (session 2, 5, 8 and 13) out of 13 flight sessions. Flight paths were reconstructed from the tracking data for each bird in the flock, with horizontal and vertical coordinates delivered by Camera 1 and coordinates in wind direction delivered by Camera 2. The data show that each bird mainly occupied a particular area in the flight section, and that this spatial preference was stable over different flight sessions. Bird Green, for example, was preferentially flying very low above the flight section’s floor, and bird Lilac preferred to fly at upwind positions in front of the flock (Fig. 1d, Extended Data Figs. 1 and 3 and Supplementary Information).Despite their preference in flight area, all birds constantly changed their spatial positions fast and rhythmically along the horizontal dimension of the flight section (Fig. 1e–g, Extended Data Figs. 2 and 4, Supplementary Video 1 and Supplementary Information). This behaviour is reminiscent of the flight behaviour of wild zebra finches: when being surprised in flight by a predator, zebra finches fly in a rapid zig-zag course low above the ground, heading for nearby vegetation16. Whether the sideways oscillating flight manoeuvres, which are performed by both wild birds in open space and domesticated birds in the wind tunnel’s flight section, are caused by the close proximity to the ground or are part of an escape reaction is yet unknown.From the tracking data, we further calculated the spatial distances in all three dimensions between all pairwise combinations of birds throughout the four flight sessions (sample rate: 24 Hz). When normalized to the maximum distance detected for each bird pairing, each dimension and each flight session, mean distances of bird pairings in all dimensions were narrowly distributed within a range of 27.7–38.0% of maximum distance (Fig. 1h and Supplementary Table 1). This may indicate that during flocking flight, zebra finches actively balance Attraction and Repulsion to maintain a stable 3D distance towards all other members of the flock. Owing to the spatial limitations in the wind tunnel’s flight section, we did not expect the zebra finches to perform large-scale flight manoeuvres with movements aligned between all flock members (Extended Data Fig. 5 and Supplementary Information), as can be observed, for example, in freely flying flocks of homing pigeons (Columba livia domestica)19 and white storks (Ciconia Ciconia)20.Visually guided horizontal repositioningWhen observing the dynamic spatial organization of our zebra finch flock, a question immediately arises: how do the birds prevent collisions during their frequent horizontal position changes? When considering the spatial limitation experienced by the flock of six birds during flight in the flight section and their highly dynamic flight style, collision rates seemed to be astonishingly low (median: 0.02 Hz; interquartile range (IQR): 0–0.03 Hz; n = 13 sessions) during flocking flight (in total 16 collisions in 13 min of analysed flight time). In birds, the visual system represents the main input channel for environmental information. To tackle the above question, we therefore first investigated the role of vision during flocking flight, and tested whether a bird’s viewing direction was correlated with the direction of horizontal position change. As gaze changes are governed by head movements in birds21, we used a bird’s head direction as an indicator for the orientation of its visual axis. We tracked (sample rate: 120 Hz) the position of a bird’s beak tip and neck in each frame of the footage during ten horizontal position changes (Fig. 2a and Supplementary Video 2) per bird, and found a strong interaction between a bird’s head angle relative to the wind direction and its direction of horizontal position change. During horizontal position changes, the birds always turned their heads in the direction of the position change (Fig. 2b). While the population’s median absolute angle of position change was 84.0° (IQR: 78.6–87.2°; n = 60) relative to 0° in wind direction, the population’s median absolute head turning angle was 36.0° (IQR: 26.4–42.5°; n = 60; see Supplementary Information for results on head movements during solo flight). The eyes of zebra finches are positioned laterally on their heads22 and each retina features a small region of highest ganglion cell density (fovea, that is, region of highest visual spatial resolution) at an area that receives visual input from horizontal positions at 60° relative to the midsagittal plane23. By turning their heads by about 36° during horizontal position changes, the zebra finches roughly align the foveal area in the retina of one eye with their direction of position change, and in the retina of the other eye with the wind direction (Fig. 2c,d). Thus, head turns in the direction of position change may indicate that the birds use visual cues while repositioning themselves within the flock. This hypothesis is supported by a study on zebra finch head movements performed during an obstacle avoidance task. In this study, instead of fixating on the obstacle, zebra finches turned their head in the direction of movement while navigating around the obstacle24.Fig. 2: Horizontal position changes are accompanied by head turns.a, Head and body orientation of bird Orange (ventral view) during one example of position changes to the right, tracked (sample rate: 120 Hz) in the footage of Camera 2. Circles: beak tip positions; plus signs: neck positions; upward pointing triangles: tail base positions. Cutouts of freeze frames of the footage taken with Camera 2 show the bird’s head and body posture for 11 time points during the position change. b, In all birds, the median angle of head turn during horizontal position change in flocking flight is positively correlated (linear mixed effects model (LMM), estimates ± s.e.m.: 2.05 ± 0.1, P  More

Neo, M. L., Eckman, W., Vicentuan, K., Teo, S.L.-M. & Todd, P. A. The ecological significance of giant clams in coral reef ecosystems. Biol. Conserv. 181, 111–123 (2015).Article 

Google Scholar 
Hill, R. W. et al. Acid secretion by the boring organ of the burrowing giant clam, Tridacna crocea. Biol. Lett. 14, 20180047 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Welsh, K., Elliot, M., Tudhope, A., Ayling, B. & Chappell, J. Giant bivalves (Tridacna gigas) as recorders of ENSO variability. Earth Planet. Sci. Lett. 307, 266–270 (2011).ADS 
CAS 
Article 

Google Scholar 
Elliot, M. et al. Profiles of trace elements and stable isotopes derived from giant long-lived Tridacna gigas bivalves: Potential applications in paleoclimate studies. Palaeogeogr. Palaeoclimatol. Palaeoecol. 280, 132–142 (2009).Article 

Google Scholar 
Killam, D., Thomas, R., Al-Najjar, T. & Clapham, M. Interspecific and intrashell stable isotope variation among the Red Sea giant clams. Geochem. Geophys. Geosyst. https://doi.org/10.1029/2019GC008669 (2020).Article 

Google Scholar 
Duprey, N., Galipaud, J.-C., Cabioch, G. & Lazareth, C. E. Isotopic records from archeological giant clams reveal a variable climate during the southwestern Pacific colonization ca. 3.0ka BP. Palaeogeogr. Palaeoclimatol. Palaeoecol. 404, 97–108 (2014).Article 

Google Scholar 
Batenburg, S. J. et al. Interannual climate variability in the Miocene: High resolution trace element and stable isotope ratios in giant clams. Palaeogeogr. Palaeoclimatol. Palaeoecol. 306, 75–81 (2011).Article 

Google Scholar 
Ayling, B. F., Chappell, J., Gagan, M. K. & McCulloch, M. T. ENSO variability during MIS 11 (424–374 ka) from Tridacna gigas at Huon Peninsula, Papua New Guinea. Earth Planet. Sci. Lett. 431, 236–246 (2015).ADS 
CAS 
Article 

Google Scholar 
Yan, H., Shao, D., Wang, Y. & Sun, L. Sr/Ca profile of long-lived Tridacna gigas bivalves from South China Sea: A new high-resolution SST proxy. Geochim. Cosmochim. Acta 112, 52–65 (2013).ADS 
CAS 
Article 

Google Scholar 
Warter, V. & Müller, W. Daily growth and tidal rhythms in Miocene and modern giant clams revealed via ultra-high resolution LA-ICPMS analysis—A novel methodological approach towards improved sclerochemistry. Palaeogeogr. Palaeoclimatol. Palaeoecol. 465, 362–375 (2017).Article 

Google Scholar 
Warter, V., Erez, J. & Müller, W. Environmental and physiological controls on daily trace element incorporation in Tridacna crocea from combined laboratory culturing and ultra-high resolution LA-ICP-MS analysis. Palaeogeogr. Palaeoclimatol. Palaeoecol. 496, 32–47 (2018).Article 

Google Scholar 
Wei, G., Sun, M., Li, X. & Nie, B. Mg/Ca, Sr/Ca and U/Ca ratios of a Porites coral from Sanya Bay, Hainan Island, South China Sea and their relationships to sea surface temperature. Palaeogeogr. Palaeoclimatol. Palaeoecol. 162, 59–74 (2000).Article 

Google Scholar 
Brahmi, C. et al. Effects of elevated temperature and pCO2 on the respiration, biomineralization and photophysiology of the giant clam Tridacna maxima. Conserv. Physiol. 9, 041 (2021).Article 
CAS 

Google Scholar 
Watson, S.-A. & Neo, M. L. Conserving threatened species during rapid environmental change: Using biological responses to inform management strategies of giant clams. Conserv. Physiol. 9, 082 (2021).
Google Scholar 
Armstrong, E. J., Dubousquet, V., Mills, S. C. & Stillman, J. H. Elevated temperature, but not acidification, reduces fertilization success in the small giant clam, Tridacna maxima. Mar. Biol. 167, 8 (2020).CAS 
Article 

Google Scholar 
Leggat, W., Buck, B. H., Grice, A. & Yellowlees, D. The impact of bleaching on the metabolic contribution of dinoflagellate symbionts to their giant clam host. Plant Cell Environ. 26, 1951–1961 (2003).CAS 
Article 

Google Scholar 
Zhou, Z., Liu, Z., Wang, L., Luo, J. & Li, H. Oxidative stress, apoptosis activation and symbiosis disruption in giant clam Tridacna crocea under high temperature. Fish Shellfish Immunol. 84, 451–457 (2019).CAS 
PubMed 
Article 

Google Scholar 
Dubousquet, V. et al. Changes in fatty acid composition in the giant clam Tridacna maxima in response to thermal stress. Biol. Open 5, 1400–1407 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Blidberg, E., Elfwing, T., Plantman, P. & Tedengren, M. Water temperature influences on physiological behaviour in three species of giant clams (Tridacnidae). In Proc. 9th International Coral Reef Symposium 561–565 (2000).Junchompoo, C., Sinrapasan, N., Penpain, C. & Patsorn, P. Changing seawater temperature effects on giant clams bleaching, Mannai Island, Rayong Province, Thailand. In Proc. Design Symposium on Conservation of Ecosystem. https://doi.org/10.13140/2.1.1906.5600 (2012).Watson, S.-A., Southgate, P. C., Miller, G. M., Moorhead, J. A. & Knauer, J. Ocean acidification and warming reduce juvenile survival of the fluted giant clam, Tridacna squamosa. Molluscan Res. 32, 177–180 (2012).
Google Scholar 
Watson, S.-A. Giant clams and rising CO2: Light may ameliorate effects of ocean acidification on a solar-powered animal. PLoS ONE 10, 1–18 (2015).CAS 

Google Scholar 
Kurihara, H. & Shikota, T. Impact of increased seawater pCO2 on the host and symbiotic algae of juvenile giant clam Tridacna crocea. Galaxea J. Coral Reef Stud. 20, 19–28 (2018).Article 

Google Scholar 
Alves Monteiro, H. J. et al. Molecular mechanisms of acclimation to long-term elevated temperature exposure in marine symbioses. Glob. Change Biol. 26, 1271–1284 (2020).ADS 
Article 

Google Scholar 
Collins, M. et al. Long-term climate change: Projections, commitments and irreversibility. In Climate Change 2013—The Physical Science Basis (eds Stocker, T. F. et al.) 1029–1136 (Cambridge University Press, 2013).
Google Scholar 
Poloczanska, E. et al. Climate change and Australian marine life. Oceanogr. Mar. Biol. 45, 407 (2007).
Google Scholar 
Ganachaud, A. S. et al. Observed and expected changes to the tropical Pacific Ocean. In Vulnerability Trop. Pac. Fish. Aquac. Clim. Change Secr. Pac. Community Noumea New Caledonia 101–187 (2011).Doney, S. C., Fabry, V. J., Feely, R. A. & Kleypas, J. A. Ocean acidification: The other CO2 problem. Annu. Rev. Mar. Sci. 1, 169–192 (2009).ADS 
Article 

Google Scholar 
Meinshausen, M. et al. The RCP greenhouse gas concentrations and their extensions from 1765 to 2300. Clim. Change 109, 213–241 (2011).ADS 
CAS 
Article 

Google Scholar 
Pierrot, D., Lewis, E. & Wallace, D. MS Excel program developed for CO2 system calculations. In ORNLCDIAC-105a Carbon Dioxide Inf. Anal. Cent. Oak Ridge Natl. Lab. US Dep. Energy Oak Ridge Tenn. Vol. 10 (2006).Mehrbach, C., Culberson, C. H., Hawley, J. E. & Pytkowicx, R. M. Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure. Limnol. Oceanogr. 18, 897–907 (1973).ADS 
CAS 
Article 

Google Scholar 
Dickson, A. G. & Millero, F. J. A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media. Deep-Sea Res. 34, 1733–1743 (1987).ADS 
CAS 
Article 

Google Scholar 
Dickson, A. G. Standard potential of the reaction: AgCl (s) + 12H2 (g) = Ag (s) + HCl (aq), and the standard acidity constant of the ion HSO4− in synthetic sea water from 273.15 to 318.15 K. J. Chem. Thermodyn. 22, 113–127 (1990).CAS 
Article 

Google Scholar 
Wolf, R. E. & Adams, M. Multi-elemental Analysis of Aqueous Geochemical Samples by Quadrupole Inductively Coupled Plasma-Mass Spectrometry (ICP-MS) 38. http://pubs.er.usgs.gov/publication/ofr20151010, https://doi.org/10.3133/ofr20151010 (2015).Schrag, D. P. Rapid analysis of high-precision Sr/Ca ratios in corals and other marine carbonates. Paleoceanography 14, 97–102 (1999).ADS 
Article 

Google Scholar 
Howell, D. C. Permutation Tests for Factorial ANOVA Designs (2009).Fox, J. & Weisberg, S. An R Companion to Applied Regression (Sage, 2019).
Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2020).Lenth, R. emmeans: Estimated Marginal Means, aka Least-Squares Means (2020).Navarro, D. Learning Statistics with R: A Tutorial for Psychology Students and other beginners (Version 0.5) (University of Adelaide, 2015).
Google Scholar 
Zhao, L., Schöne, B. R. & Mertz-Kraus, R. Controls on strontium and barium incorporation into freshwater bivalve shells (Corbicula fluminea). Palaeogeogr. Palaeoclimatol. Palaeoecol. 465, 386–394 (2017).Article 

Google Scholar 
Bragg, W. L. The structure of some crystals as indicated by their diffraction of X-rays. Proc. R. Soc. Lond. Ser. Contain. Pap. Math. Phys. Character 89, 248–277 (1913).ADS 
CAS 

Google Scholar 
Bragg, W. L. The structure of aragonite. Proc. R. Soc. Lond. Ser. Contain. Pap. Math. Phys. Character 105, 16–39 (1924).ADS 
CAS 

Google Scholar 
Killam, D., Al-Najjar, T. & Clapham, M. Giant clam growth in the Gulf of Aqaba is accelerated compared to fossil populations. Proc. R. Soc. B Biol. Sci. 288, 20210991 (2021).CAS 
Article 

Google Scholar 
Waters, C. G. Biological Responses of Juvenile Tridacna maxima (Mollusca: Bivalvia) to Increased pCO2 and Ocean Acidification (The Evergreen State College, 2008).
Google Scholar 
Toonen, R. J., Nakayama, T., Ogawa, T., Rossiter, A. & Delbeek, J. C. Growth of cultured giant clams (Tridacna spp.) in low pH, high-nutrient seawater: Species-specific effects of substrate and supplemental feeding under acidification. J. Mar. Biol. Assoc. U. K. 92, 731–740 (2012).CAS 
Article 

Google Scholar 
Hart, A. M., Bell, J. D. & Foyle, T. P. Growth and survival of the giant clams, Tridacna derasa, T. maxima and T. crocea, at village farms in the Solomon Islands. Aquaculture 165, 203–220 (1998).Article 

Google Scholar 
Van Wynsberge, S. et al. Growth, survival and reproduction of the giant clam Tridacna maxima (Röding 1798, Bivalvia) in two contrasting lagoons in French Polynesia. PLoS ONE 12, 1–20 (2017).
Google Scholar 
Lucas, J. S., Nash, W. J., Crawford, C. M. & Braley, R. D. Environmental influences on growth and survival during the ocean-nursery rearing of giant clams, Tridacna gigas (L.). Aquaculture 80, 45–61 (1989).Article 

Google Scholar 
Schwartzmann, C. et al. In situ giant clam growth rate behavior in relation to temperature: A one-year coupled study of high-frequency noninvasive valvometry and sclerochronology. Limnol. Oceanogr. 56, 1940–1951 (2011).ADS 
Article 

Google Scholar 
Syazili, A., Syafiuddin, N. A. & Jompa, J. Effect of ocean acidification and temperature on growth, survival, and shell performance of fluted giant clams (Tridacna squamosa). IOP Conf. Ser. Earth Environ. Sci. 473, 012141 (2020).Article 

Google Scholar 
Li, J. et al. Assessment of the juvenile vulnerability of symbiont-bearing giant clams to ocean acidification. Sci. Total Environ. 812, 152265 (2022).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Li, S. et al. Cloning and expression of a pivotal calcium metabolism regulator: Calmodulin involved in shell formation from pearl oyster (Pinctada fucata). Comp. Biochem. Physiol. B Biochem. Mol. Biol. 138, 235–243 (2004).PubMed 
Article 
CAS 

Google Scholar 
Wang, X., Li, C., Lv, Z., Zhang, Z. & Qiu, L. A calcification-related calmodulin-like protein in the oyster Crassostrea gigas mediates the enhanced calcium deposition induced by CO2 exposure. Sci. Total Environ. 833, 155114 (2022).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Rees, T., Fitt, W. & Yellowlees, D. The haemolymph and its temporal relationship with zooxanthellae metabolism in the giant clam symbiosis [Conference paper]. In ACIAR Proc.-Aust. Cent. Int. Agric. Res. Aust. (1993).Leggat, W., Rees, T. A. V. & Yellowlees, D. Meeting the photosynthetic demand for inorganic carbon in an alga-invertebrate association: Preferential use of CO2 by symbionts in the giant clam Tridacna gigas. Proc. Biol. Sci. 267, 523–529 (2000).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ip, Y. K. et al. Molecular characterization, light-dependent expression, and cellular localization of a host vacuolar-type H+-ATPase (VHA) subunit A in the giant clam, Tridacna squamosa, indicate the involvement of the host VHA in the uptake of inorganic carbon and. Gene 659, 137–148 (2018).CAS 
PubMed 
Article 

Google Scholar 
Armstrong, E. J., Roa, J. N., Stillman, J. H. & Tresguerres, M. Symbiont photosynthesis in giant clams is promoted by V-type H+-ATPase from host cells. J. Exp. Biol. https://doi.org/10.1242/jeb.177220 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Sano, Y. et al. Past daily light cycle recorded in the strontium/calcium ratios of giant clam shells. Nat. Commun. 3, 761 (2012).ADS 
PubMed 
Article 
CAS 

Google Scholar 
Adams, A. L., Needham, E. W. & Knauer, J. The effect of shade on water quality parameters and survival and growth of juvenile fluted giant clams, Tridacna squamosa, cultured in a land-based growth trial. Aquac. Int. 21, 1311–1324 (2013).CAS 
Article 

Google Scholar 
Rossbach, S., Saderne, V., Anton, A. & Duarte, C. M. Light-dependent calcification in Red Sea giant clam Tridacna maxima. Biogeosciences 16, 2635–2650 (2019).ADS 
CAS 
Article 

Google Scholar 
Ip, Y. K. et al. The whitish inner mantle of the giant clam, Tridacna squamosa, expresses an apical plasma membrane Ca2+-ATPase (PMCA) which displays light-dependent gene and protein expressions. Front. Physiol. 8, 781 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Berner, R. A. The role of magnesium in the crystal growth of calcite and aragonite from sea water. Geochim. Cosmochim. Acta 39, 489–504 (1975).ADS 
CAS 
Article 

Google Scholar 
Alibert, C. et al. Source of trace element variability in Great Barrier Reef corals affected by the Burdekin flood plumes. Geochim. Cosmochim. Acta 67, 231–246 (2003).ADS 
CAS 
Article 

Google Scholar 
McCulloch, M. et al. Coral record of increased sediment flux to the inner Great Barrier Reef since European settlement. Nature 421, 727–730 (2003).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Sinclair, D. & Mcculloch, M. Corals record low mobile barium concentrations in the Burdekin River during the 1974 flood: Evidence for limited Ba supply to rivers? Palaeogeogr. Palaeoclimatol. Palaeoecol. 214, 155–174 (2004).Article 

Google Scholar 
Fleitmann, D. et al. East African soil erosion recorded in a 300 year old coral colony from Kenya. Geophys. Res. Lett. 34, L04401 (2007).ADS 
Article 

Google Scholar 
Prouty, N. G., Field, M. E., Stock, J. D., Jupiter, S. D. & McCulloch, M. Coral Ba/Ca records of sediment input to the fringing reef of the southshore of Moloka’i, Hawai’i over the last several decades. Mar. Pollut. Bull. 60, 1822–1835 (2010).CAS 
PubMed 
Article 

Google Scholar 
Fallon, S. J., McCulloch, M. T., van Woesik, R. & Sinclair, D. J. Corals at their latitudinal limits: Laser ablation trace element systematics in Porites from Shirigai Bay, Japan. Earth Planet. Sci. Lett. 172, 221–238 (1999).ADS 
CAS 
Article 

Google Scholar 
Reuer, M. K., Boyle, E. A. & Cole, J. E. A mid-twentieth century reduction in tropical upwelling inferred from coralline trace element proxies. Earth Planet. Sci. Lett. 210, 437–452 (2003).ADS 
CAS 
Article 

Google Scholar 
Montaggioni, L. F., Le Cornec, F., Corrège, T. & Cabioch, G. Coral barium/calcium record of mid-Holocene upwelling activity in New Caledonia, South-West Pacific. Palaeogeogr. Palaeoclimatol. Palaeoecol. 237, 436–455 (2006).Article 

Google Scholar 
Ourbak, T. et al. A high-resolution investigation of temperature, salinity, and upwelling activity proxies in corals: Activity proxies in corals. Geochem. Geophys. Geosyst. 7, 1. https://doi.org/10.1029/2005GC001064 (2006).CAS 
Article 

Google Scholar 
Alibert, C. & Kinsley, L. A 170-year Sr/Ca and Ba/Ca coral record from the western Pacific warm pool: 1. What can we learn from an unusual coral record? J. Geophys. Res. 113, C04008 (2008).ADS 

Google Scholar 
Alibert, C. & Kinsley, L. A 170-year Sr/Ca and Ba/Ca coral record from the western Pacific warm pool: 2. A window into variability of the new ireland coastal undercurrent. J. Geophys. Res. 113, C06006 (2008).ADS 

Google Scholar 
Agbaje, O. B. A. et al. Architecture of crossed-lamellar bivalve shells: The southern giant clam (Tridacna derasa, Röding, 1798). R. Soc. Open Sci. 4, 170622 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Clark, M. S. et al. Deciphering mollusc shell production: The roles of genetic mechanisms through to ecology, aquaculture and biomimetics. Biol. Rev. https://doi.org/10.1111/brv.12640 (2020).Article 
PubMed 

Google Scholar 
Wilkerson, F. P. & Trench, R. K. Uptake of dissolved inorganic nitrogen by the symbiotic clam Tridacna gigas and the coral Acropora sp.. Mar. Biol. 93, 237–246 (1986).CAS 
Article 

Google Scholar 
Summons, R. E., Boag, T. S. & Osmond, C. B. The effect of ammonium on photosynthesis and the pathway of ammonium assimilation in Gymnodinium microadriaticum in vitro and in symbiosis with tridacnid clams and corals. Proc. R. Soc. Lond. B Biol. Sci. 227, 147–159 (1986).ADS 
CAS 
Article 

Google Scholar 
Onate, J. & Naguit, M. A preliminary study on the effect of increased nitrate concentration on the growth of giant clams Hippopus hippopus. In Cult. Giant Clams Bivalvia Tridacnidae Aust. Cent. Int. Agric. Res. Canberra 57–61 (1989).Hastie, L. C., Watson, T. C., Isamu, T. & Heslinga, G. A. Effect of nutrient enrichment on Tridacna derasa seed: Dissolved inorganic nitrogen increases growth rate. Aquaculture 106, 41–49 (1992).CAS 
Article 

Google Scholar 
Belda, C. A., Lucas, J. S. & Yellowlees, D. Nutrient limitation in the giant clam-zooxanthellae symbiosis: Effects of nutrient supplements on growth of the symbiotic partners. Mar. Biol. 117, 655–664 (1993).Article 

Google Scholar 
Belda-Baillie, C., Leggat, W. & Yellowlees, D. Growth and metabolic responses of the giant clam-zooxanthellae symbiosis in a reef-fertilisation experiment. Mar. Ecol. Prog. Ser. 170, 131–141 (1998).ADS 
CAS 
Article 

Google Scholar 
Calosi, P. et al. Multiple physiological responses to multiple environmental challenges: An individual approach. Integr. Comp. Biol. 53, 660–670 (2013).CAS 
PubMed 
Article 

Google Scholar 
Tanner, R. L. & Dowd, W. W. Inter-individual physiological variation in responses to environmental variation and environmental change: Integrating across traits and time. Comp. Biochem. Physiol. A. Mol. Integr. Physiol. 238, 110577 (2019).CAS 
PubMed 
Article 

Google Scholar 
Guscelli, E., Spicer, J. I. & Calosi, P. The importance of inter-individual variation in predicting species’ responses to global change drivers. Ecol. Evol. 9, 4327–4339 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Telesca, L. et al. Biomineralization plasticity and environmental heterogeneity predict geographical resilience patterns of foundation species to future change. Glob. Change Biol. 25, 4179–4193 (2019).ADS 
Article 

Google Scholar 
Yan, H., Shao, D., Wang, Y. & Sun, L. Sr/Ca differences within and among three Tridacnidae species from the South China Sea: Implication for paleoclimate reconstruction. Chem. Geol. 390, 22–31 (2014).ADS 
CAS 
Article 

Google Scholar 
Warter, V., Mueller, W., Wesselingh, F. P., Todd, J. A. & Renema, W. Late Miocene seasonal to subdecadal climate variability in the Indo-West Pacific (East Kalimantan, Indonesia) preserved in giant clams. Palaios 30, 66–82 (2015).ADS 
Article 

Google Scholar 
Gannon, M. E., Pérez-Huerta, A., Aharon, P. & Street, S. C. A biomineralization study of the Indo-Pacific giant clam Tridacna gigas. Coral Reefs 36, 503–517 (2017).ADS 
Article 

Google Scholar 
Zhao, L. et al. A review of transgenerational effects of ocean acidification on marine bivalves and their implications for sclerochronology. Estuar. Coast. Shelf Sci. 235, 106620 (2020).CAS 
Article 

Google Scholar  More

Sample collection from sea turtlesIn total, 45 serum samples from 33 juvenile green turtles (C. mydas), including 6 sea turtles with tumors, 5 juvenile hawksbill turtles (Eretmochelys imbricate), and 7 olive ridley turtles (Lepidochelys olivacea) (juvenile = 5; sub-adult = 2). All turtles were sourced from: eastern Taiwan (n = 24), southern Taiwan (n = 14), central Taiwan (n = 6), and northern Taiwan (n = 1). Among the 45 sea turtle samples, 6 green turtles developed FP (n = 1 with tumor score 1; n = 1 with tumor score 2; n = 4 with tumor score 3)32, while 39 did not have FP. FP tumor tissues were collected from 6 green turtles (from shoulder/flippers/inguinal regions) with FP during surgical procedures. Regarding the collection of normal skin tissues, one normal skin tissue (from shoulder) was collected from one necropsied dead green turtles (stranding and discovered from southern Taiwan) confirmed without FP. All tissue samples were fixed in 10% neutral buffered formalin prior to further analysis. In this study, all sea turtles were discovered and rescued through the official reporting system of the Marine Animal Rescue Network (established by the Ocean Conservation Administration) and admitted to the National Museum of Marine Biology and Aquarium (NMMBA), between 2017 and 2020.Detection of ChHV5 DNA by polymerase chain reaction (PCR)Total DNA was extracted from blood of 45 sea turtles by DNeasy blood & tissue kit (Cat. No. 69504, Qiagen, Valencia, CA, USA) following manufacturer’s instructions. Subsequently, the ChHV5 infection status all 45 sea turtles was determined by PCR using primers targeting on UL18 (capsid protein gene), UL22 (glycoprotein H gene), and UL27 (glycoprotein B gene) regions4. The sequence of primer sets are: UL18F: 5′-CACCACGAGGGGGAAAATGA, UL18R:5′-TCAAATCCCCCGTTCACTCG; UL22F: 5′-ACGGCGTTGGCTAGTGAATC, UL22R: 5′-GCAGTTCGGTACACACCTCT; UL27F: 5′-TAACAAGAAAGAACCGCGCG; UL27R: 5′-ATTTTCCCGGTCAGTGCCAA. PCR amplifications were performed in a total volume of 50 μl. The reaction included 1 μl of the template DNA, 1 μl of each primer (10 μM), 22 μl of distilled water (DDW), and 25 μl of the AmpliTaq Gold® 360 Master Mix (Cat. No. 4398876, Life Technologies, Valencia, CA, USA). The thermocycle for amplification was: Initial denaturing at 95 °C for 10 min, followed by 40 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 60 s, and then a final extension at 72 °C for 7 min. Results were visualized by gel electrophoresis (2% agarose) with SYBR Safe DNA Gel Stain (Cat. No. S33102, Invitrogen, Carlsbad, CA, USA).Sequence optimization of the UL27 gene for expression of the ChHV5 glycoprotein protein using E. coli
To express large quantities of ChHV5 gB, we adopted the prokaryotic Escherichia coli (E. coli) expression system. The construct (namely UL27/pUC57) containing sequences of the full length UL27 fused with FLAG tag sequence (GenBank accession no. AF035003.3) was synthesized by Allbio Science Co., Ltd, Taiwan. The sequence information of the glycoprotein (gB) datasets used and analyzed for protein expression during the current study was obtained and available from the GenBank repository [https://www.ncbi.nlm.nih.gov/nuccore/AF035003.3]. Considering the difference in tRNA-codon usages between prokaryotes and eukaryotes would possibly affect subsequent protein expression, the optimized UL27 gene sequence, without altering the translated amino acid sequences, to fit the E. coli expression system was synthesized. The codon optimized UL27 gene was further used as the template for amplification of different gene fragments by Polymerase Chain Reaction (PCR).Construction of plasmids expressing partial fragments of ChHV5 gB proteinTo determine the relative antigenicity and also to increase the expression yield, plasmids expressing various regions of gB protein were constructed. Briefly, the five regions covering different fragments of the UL27 gene were amplified from plasmid UL27/pUC57 by PCR using specific primer sets with built-in restriction enzyme sequences shown as underlined in Table 1. The thermal cycling conditions were: 98 °C (5 min) followed by 35 cycles of denaturation (98 °C, 30 s), annealing (58 °C, 1 min), and extension (72 °C, 2 min), and finished with a final extension (72 °C, 10 min). PCR amplicons with expected sizes were isolated from gel and trimmed with the restriction enzymes followed by ligation with vectors either pET24a or pET32b (Novagen, Germany) linearized with the same restriction enzymes. The resulting plasmids with expected insert sizes as confirmed by restriction enzymes were sent for automated DNA sequencing (Mission Biotech, Taipei, Taiwan).Table 1 Information on the constructs expressing the UL27 fragments. The bold characters indicate sequences recognized by restriction enzymes for the ease of further cloning procedure.Full size tableExpression of recombinant gB fragments in E. coli
In the current study, the recombinant gB protein is a key reagent that served as antigen for seroprevalence of ChHV5 as well as for the generation of ChHV5 gB antibody (conducted by Yao-Hong Biotech Inc., Taiwan). The plasmids expressing individual gB fragment were transformed into E. coli host cells, strain BL21 (DE3), Rosetta. Expression of all the recombinant gB fragments was induced by 0.8 mM of isopropyl β-d-1-thiogalactopyranoside (IPTG) at 28 °C for 16 h. As all the gB fragments cloned into the pET series vectors were expressed as a fusion protein with a 6-histidine tag at C-terminus end, they could be further purified by Ni–NTA column chromatography using the chelating Sepharose Fast Flow (GE Healthcare) following the method described in one previous study33. The yield and purity of recombinant gB proteins were confirmed by sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE). Subsequently, 6 M urea and 0.4 M imidazole contained in the purified protein were depleted by step-wise dialysis against 1 × PBS buffer (0.02 M phosphate, 0.15 M NaCl) with gradually decreased concentrations of urea at 4 °C. The concentration of recombinant proteins were then estimated by National Institutes of Health ImageJ software (https://imagej.nih.gov/ij/, 1997–2018.) using the standard curve established by bovine serum albumin (BSA) with known concentrations42.Western blot analysisRecombinant gB fragments were separated by 12.5% or 15% SDS-PAGE and electrotransferred to PVDF membrane by using Mini Proten III apparatus (Cat. No. 165-3301, BioRad). The filters were blocked in PBS-T buffer (0.02 M phosphate, 0.15 M NaCl, 0.05% Tween-20) containing 5% skimmed milk and reacted with mouse anti-his tag antibody (1:5,000, Cat. No. GTX40628, GeneTex) at 4 °C for overnight. After six-time wash with PBS-T buffer, the PVDF filter was then incubated with the secondary antibody, 1:5000 diluted goat anti-mouse IgG conjugated with horseradish peroxidase (HRP), or 1:500 diluted Protein A/G-HRP (Cat. No. 32400, Thermo fisher scientific™, United States) for sea turtle antibody detection, at room temperature for 1 h followed by PBS-T wash to remove the unbound antibodies. Ultimately, the signal was detected by ECL reagents (Thermo Fisher Scientific, United States) and the image was acquired by ImageQuant LAS 4000 Mini (GE Healthcare).Immunohistochemical (IHC) analysisTo establish IHC protocol, normal skin tissue from PCR-negative sea turtles served as the negative control. In total, the FP on skin tissue from six individual sea turtles that were detected positive for ChHV5 DNA (positive tissue samples), and one normal tissue detected negative (the negative tissue) were included in the IHC analysis.IHC procedure was conducted as reported in our previous study34. In brief, sections of formalin-fixed and wax-embedded skin tissues of sea turtles were made using a rotary microtome (Leica RM2245, Leica Biosystems, Germany) and were further deparaffinized and rehydrated. Antigen retrieval was carried out by heat-induced epitope retrieval method: slides immersed into boiled sodium citrate buffer (0.01 M, pH 6.0), which was preheated up to 100 °C, for 20 min and cooled at room temperature for 20 min. Subsequently, the slides were incubated with peroxidase-blocking reagent (Cat. No. S200389, Dako, Denmark) for 30 min, and then treated with or without primary antibodies (the anti-gB serum prepared from this study). In each interval of the following procedures, sections were rinsed with a mixture of TBST buffer. Tissue sections were then reacted with secondary antibody (HRP anti-rabbit/mouse, DAKO, Denmark), followed by incubation of DAB and chromogen (dilution 1 μL in 100 μL) from a commercial ChemMate EnVision detection kit (Cat. No. K5007, Dako, Denmark). Ultimately, tissue sections were counterstained with Mayer’s hematoxylin reagents (Code S3309, Dako, Denmark) for 2 min followed by wash with DDW, and reacted with 37 mM ammonia water for 5 s and rinsed with DDW.Immunofluorescent assay (IFA)Human 293 T cells were transfected with plasmids expressing full-length ChHV5 gB protein fused with FLAG tag at its C-terminus. At 24 h post transfection, 293 T cells (CRL-3216, ATCC, USA) were fixed with 2% formaldehyde for 10 min, followed by permeabilization with 0.1% Triton X-100 for another 10 min. Subsequently, cells were incubated with anti-FLAG antibody (1:500) (F7425; Sigma-Aldrich), or antisera (F1, F2, F3, F2–3) at the dilution of 1:500 for 1 h at room temperature. After six times of washes with PBS containing1% bovine serum, goat anti-mouse IgG (1:2,000 fold diluted) (Cat. No. A28175, Alexa Fluor® 488, Invitrogen) was used as secondary antibody. After one-hour incubation, nuclei were stained with 4, 6-diamidino-2-phenylindole (DAPI, Cat. No. D9542, Sigma-Aldrich) for 10 min, followed by confocal microscopy (FV1000, Olympus, Tokyo, Japan) with Olympus FV10-ASW 1.3 viewer software.Statistical analysisTo evaluate the association between seropositivity and FP or viremia tested by PCR of UL27 gene, Fisher’s exact test was applied due to very limited number of sea turtles with FP. The statistical significance was determined by p  More

Prosopis juliflora leaves collection and processing for RibotypingProsopis juliflora species of the genus Prosopis, family of Fabaceae had its genetic variation in Doha evaluated. Seven samples of P. juliflora leaves were collected from six different locations in Doha, Qatar, during five field trips. Plant leaves were collected after proper permissions and all methods were carried out in accordance with relevant guidelines and regulations. Trees in all locations were naturally growing around urbanization areas in their normal arid habitat without artificial irrigation, samples were collected from fully mature trees. Table 1 shows the samples details. Figure 1 shows the location sites of where the samples were collected on the map of Qatar, Doha. Leaf samples were kept in sterile labeled bags until having reached the laboratory where few leaflets were washed with sterile distilled water and sterilized using 70% ethanol to be used for DNA extraction.Table 1 Location details of the collection sites of P. juliflora leaves.Full size tableFigure 1Location map of collection sites of P. juliflora leaf samples (ArcGIS software).Full size imageRibotyping analysisThe leaflets of each sample were transferred into a sterile mortar previously cooled at -20 ˚C and used for DNA extraction following the kit manufacturer instructions (DNeasy Plant Mini Kit-QIAGEN-USA).Extracted DNA of each sample were subject to PCR using ITS1 and ITS4 primers. PCR products obtained were purified using the Invitrogen Quick PCR Purification Kit (QIAGEN, Germany) as indicated by the manufacturer and sequenced using Sanger sequencer (3130/3130xl DNA Analyzers, Thermofisher Scientific, USA) as previously described22.Sanger sequencer raw data was read using BioEdit software. Basic Local Alignment Search Tool (BLAST) network services of the National Centre for Biotechnology Information (NCBI) database were used to compare the obtained sequences to the existing sequences. Sequences were submitted to NCBI for accession numbers. The various P. juliflora ribosomal sequences obtained were also uploaded on MEGA-X software and the phylogeny tree was generated using the neighbor-joining algorithm26.Minimum inhibitory concentrations of PJ-WS-LE extracts prepared using leaf samples collected from various locations against A. alternata and C. gloeosporioides
Preparation of PJ-WS-LE extractFresh, young full leaves of P. juliflora were collected from various locations as indicated in Fig. 1. Samples were washed, dried and ground into powder to be used in the preparation of PJ-WS-LE extract as previously described22. Briefly, every 20 g of the leaf powder were incubated in 200 mL of 70% ethanol for 48 h. The supernatant has its solvent evaporated, the extract was then re-dissolved in sterile distilled water. Only water-soluble phytochemicals were tested by centrifuging the final preparation tubes and excluding the pellet. Stock solution of 100 g L−1 was stored at 4 °C to be used for later experiments. PJ-WS-LE extract concentration used in treatments was 8 g L−1 which is double the highest minimum inhibitory concentration of the extract against spoiling microorganisms as previously determined22.Determination of minimum inhibitory concentrationThe MIC test was conducted in a sterile 96-well plate, with each well containing 100 μl of potato dextrose broth (PDB) (HIMEDIA-India). Every four wells made one replication, nine different concentrations of the crude extracts were tested (1:1 dilutions) ranging from 42 to 0.16 g L−1. Wells were then inoculated with one of the two tested microorganisms’ spore suspensions (A. alternata and C. gloeosporioides). The last three rows are control rows: no spores and no extract control wells, negative control with spores but no extract wells, and positive control with spores and 10 µl of the fungicidal Clatrimazole (1%) wells.Fungal spore suspensions were adjusted to the range of 104 spores L−1 using a 10 day old fungal plate and sterile distilled water, the spore concentration was calculated using a heamatocytometer.Fungal growth in each well was monitored using Resazurin (HIMEDIA-India) dye. Upon cells division, Resazurin changes its color from blue to pink and fluorescent27. Results were taken within 48 h of incubation at 25 °C. MIC was recorded as the last extract concentration that shows no change in the color of Resazurin within the incubation period.Curative and preventive effects of PJ-WS-LE extract against A. alternata and C. gloeosporioides induced infection in mangoesPathogensThe two fungal strains used C. gloeosporioides and A. alternata were obtained from our laboratory collection, Department of Biological and Environmental Sciences, Qatar University, Qatar. Both fungal isolates were previously isolated from locally collected fruit samples. Isolates were molecularly identified by sequencing the Internal Transcribed Spacer (ITS) regions of fungal ribosomal DNA (rDNA) that was amplified by PCR. Identified fungal isolates were given the strains code of AaltQU17 for A. alternata and CgloQU17 for C. gloeosporioides22. Preserved cultures were sub-cultured on potato dextrose agar (PDA) plates and incubated at 25 °C for 10 days. Plates were then flooded with 10 mL of sterile distilled water each, to prepare the needed spores suspension solutions. The concentrations of spores suspensions were adjusted to 106 spores L−1 using a heamatocytometer18.FruitThe mango (Mangifera indica) type known as Neelam imported from India was used in the experiments. Fruit were bought from the whole sale market upon their arrival to the country. Only undamaged mature fruit were used in the experiment. Fruits chosen were ripen but not yet soft with firmness average of 20 ± 5.1 N, weight average of 177.61 ± 0.2 g and TSS average of 70 ± 5.3%. Fruit were first washed with sink water and sterilized twice with 70% ethanol to be then washed with sterile distilled water and left to air dry.Preventive and curative effects of PJ-WS-LE extractWounded mango fruit were used during the experiment, the wounds were made through three needle pricks (2 mm deep) in three different places for each plant using a sterile syringe. A completely randomized design was used and each treatment was made of a triplicate of 10 fruit each. The experiment was repeated twice.PJ-WS-LE extract of leaves collected from Qatar university field was first tested for its efficacy in preventing fungal contamination in wounded mango fruit (preventive effect). Therefore, the wounded zone of each fruit was sprayed with 8 g L−1 PJ-WS-LE extract and then left to air-dry. Once dry the fruit were sprayed again with the extract at the same concentration and left to dry. Control fruit were only treated with sterile distilled water without the plants extract. After two hours all wounds were inoculated with 20 μL of conidia aqueous solution (106 spores mL−1) of one of the tested fungi. The extract was then tested for its ability to cure fungal contamination in wounded fruit. Therefore, wounds were inoculated first with 20 μL of conidia aqueous solution (106 spores mL−1) and left to dry. Wounds were then sprayed twice with 8 g L−1 PJ-WS-LE extract.All mangoes were stored in sterilized plastic trays inside an incubator at 25 °C and 75% humidity. Fruit were observed every 24 h for 5 days for C. gloeosporioides inoculated fruit and for 10 days for A. alternata inoculated fruit. Three parameters were recorded at the end of the experiment: disease incidence (DI), disease severity (DS), and percent plant extract efficacy (%EE). To calculate disease severity, the diameter of the infected area of each fruit was measured in two perpendicular directions and mean diameter mycelial growth was calculated28,29.$$mathrm{DI}=frac{(mathrm{Number, of, rotten, fruit})times 100}{mathrm{Total, number, of, fruit}}$$$$mathrm{DS }=frac{(mathrm{Average, lesion, diameter, of, treated, plants})times 100}{mathrm{Average, Lesion, diameter, of, control, plants})}$$$$mathrm{%EE}=frac{(mathrm{Disease, incidence, in, Control, batch}-mathrm{Disease, incidence, in, treated, batch})times 100}{mathrm{Disease, incidence, in, Control, batch}}$$End of the trial samples firmnessAt the end of the trial, remaining mango fruit were tested for their flesh quality using a penetrometer (Agriculture Solutions, USA) to test the flesh firmness. Fruit were peeled, then the stainless steel probe of the instrument was inserted in three different points towards the equator of the fruit. Firmness in Newton was recorded and compared with standard fruit firmness to judge fruit quality18.Effectiveness of PJ-WS-LE extract as long-term coating material and the preservative value of its chitosan-embedded formCoating solutions preparationChitosan solution of 1% concentration was prepared by stirring chitosan powder (CAS 9012-76-4, Himedia, India) in 1% glacial acetic acid (IsoLab, Germany) overnight. The final chitosan solution pH was adjusted to 5.6 using 0.1 M NAOH (Sigma-Aldrich, Germany). To prepare PJ-WS-LE extract chitosan-embedded coating material, filter-sterilized PJ-WS-LE extract stock solution was added to 1% chitosan to achieve a final concentration of 8 g L−130.Samples preparationEighty-four mango samples chosen as described above, were divided into four groups of 18 samples each. Samples were divided into four treatment batches and treated as following:

Batch A: non-treated fruit.

Batch B: PJ-WS-LE extract at 8 g L–1 was used to spray the fruit.

Batch C: 1% chitosan was used to spray the fruit.

Batch D: 8 g L−1 PJ-WS-LE extract embedded in 1% chitosan was used to spray the fruit.

Every experimental replicate was made up of three mango samples that were stored together in one sterile bag at 4 °C. The number of replications per treatment was six. The experiment was repeated twice31.Evaluation of sensory qualityA five-points scale was used for the evaluation of the sensory quality of the samples for overall quality, smell, and color change. The attributes were evaluated weekly using the fruit of one experimental replicate. Scores were given using the following scale: 5 points indicate “extremely liked”, 4 points indicate “liked”, 3 points indicate “acceptable” 2 points indicate “disliked” and 1 point indicates “extremely disliked”. The weekly average score per batch was also calculated32.Estimation of weight lossUpon treatment at day zero, all mango samples were weighed and their weights were recorded as initial weights. Weights of all remaining samples were measured at the end of every week. The variation between the start weight and weekly weights is calculated as weekly weight loss. The average percent of weekly weight loss of each batch was calculated32.Determination of samples firmnessThe samples of each experimental replicate evaluated on a weekly basis had their firmness measured as previously described. The weekly average samples firmness (N) of every treatment batch was also calculated33.pH measurementMango fruit of each experimental replicate were blended weekly into juice, after filtration, a digital pH meter (Jenway, UK) was used to measure pH. The weekly average fruit pH of every treatment batch was also calculated. The pH meter was calibrated using a buffer solution of pH 734.Total soluble solids (TSS) measurementTotal soluble solids of the prepared mango juice samples were measured in percent brix using a refractometer (ANTAHI, New Zealand). The weekly average fruit TSS (%) for each treatment batch was also calculated. The refractometers was calibrated using distilled water35.DPPH radical scavenging assayA 1/10 mango juice dilution was prepared using sterile distilled water. 100 μL of each dilution was mixed with 1 mL of 2,2-diphenyl-1-picrylhydrazyl (DPPH) (100 mg L−1) to be incubated in the dark at 37 °C for 45 min. After incubation, samples were centrifuged and the pellet was discarded. The intensity of the change in color of the supernatant was measured by spectrophotometry at 517 nm using methanol as a blank. 100 μL of methanol in 1 mL DPPH was used as the control for the experiment. Percent radical scavenging activity was calculated as per the below formula:$$ % {text{ radical scavenging activity}}, = ,left( {{text{absorbance of the control solution}} – {text{ absorbance of the juice sample}}} right)*{1}00/{text{absorbance of the control solution}}. $$The weekly average % radical scavenging activity for each treatment batch was finally calculated31.Statistical analysisThe experimental design used was Completely Randomized Design (CRD). One-way ANOVA followed by Tukey Post-Hoc test was used to evaluate the significance of the weekly percent change in weight among treatment batches at P ≤ 0.05. The significances of pH and TSS variation within different treatment batches were evaluated using One-way ANOVA test at P ≤ 0.05. Data was presented as average ± standard error of the Means (SEM). SPSS (Ver. 27, SPSS Inc. Chicago, USA) was used to perform the statistical analysis tests. More

Printzen, C. & Lumbsch, H. T. Molecular evidence for the diversification of extant lichens in the late cretaceous and tertiary. Mol. Phylogenet. Evol. 17, 379–387 (2000).CAS 
PubMed 
Article 

Google Scholar 
Kraichak, E. et al. A Tale of two Hyper-diversities: Diversification dynamics of the two largest families of lichenized fungi. Sci. Rep. https://doi.org/10.1038/srep100288 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Leavitt, S. D., Lumbsch, H. T., Stenroos, S. & Clair, L. L. S. Pleistocene speciation in North American lichenized fungi and the the impact of alternative species circumscriptions and rates of molecular evolution on divergence estimates. PLoS ONE 8, e85240 (2013).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Gaya, E. et al. The adaptive radiation of lichen-forming Teloschistaceae is associated with sunscreening pigments and bark-to-rock substrate shift. PNAS 112, 11600–11605 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Schneider, K., Resl, P. & Spribille, T. Escape from the cryptic species trap: Lichen evolution on both sides of a cyanobacterial acquisition event. Mol. Ecol. 25, 3453–3468 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Widhelm, T. J. et al. Oligocene origin and drivers for diversification in the genus Sticta (Lobariaceae, Ascomycota). Mol. Phylogenetic Evol. 126, 58–73 (2018).Article 

Google Scholar 
Vamosi, J. C. & Vamosi, S. M. Factors influencing diversification in angiosperms: At the crossroads of intrinsic and extrinsic traits. Am. J. Bot. 98, 460–471 (2011).PubMed 
Article 

Google Scholar 
Wagner, C. E., Harmon, L. J. & Seehausen, O. Ecological opportunity and sexual selection together predict adaptive radiation. Nature 487, 366–369 (2012).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Karunarathne, P. et al. Intraspecific ecological niche divergence and reproductive shifts foster cytotype displacement and provide ecological opportunity to polyploids. Ann. Bot. 121, 1183–1196 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Nakov, T., Beaulieu, J. & Alverson, A. Accelerated diversification is related to life history and locomotion in a hyperdiverse lineage of microbial eukaryotes (Diatoms, Bacillariophyta). New Phytol. 219, 462–473 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Tripp, E. A. Is asexual reproduction an evolutionary dead end in lichens?. Lichenologist 48, 559–580 (2016).Article 

Google Scholar 
Tripp, E. A. & Lendemer, J. C. Twenty-seven modes of reproduction in the obligate lichen symbiosis. Brittonia 70, 1–14 (2018).Article 

Google Scholar 
Bowler, P. A. & Rundell, P. W. Reproductive strategies in lichens. Bot. J. Linn. Soc. 70, 325–340 (1975).Article 

Google Scholar 
Honegger, R. Developmental biology of lichens. New Phytol. 125, 659–677 (1993).PubMed 
Article 

Google Scholar 
Buschbom, J. & Mueller, G. M. Testing “Species Pair” hypotheses: Evolutionary processes in the lichen-forming species complex Porpidia flavocoerulescens and Porpidia melinodes. Mol. Biol. Evol. 23, 574–586. https://doi.org/10.1093/molbev/msj063 (2006).CAS 
Article 
PubMed 

Google Scholar 
Sanders, W. B. Complete life cycle of the lichen fungus Calopadia puiggarii (Pilocarpaceae, Ascomycetes) documented in situ: Propagule dispersal, establishment of symbiosis, Thallus development, and formation of sexual and asexual reproductive structures. Am. J. Bot. 101, 1836–1848 (2014).PubMed 
Article 

Google Scholar 
Poelt, J. Flechtenflora und eiszeit in Europa. Phyton (Horn) 10, 206–214 (1963).
Google Scholar 
Stofer, S. et al. Species richness of lichen functional groups in relation to land use intensity. Lichenologist 38, 331–353 (2006).Article 

Google Scholar 
Ludwig, L. R., Summerfield, T. C., Lord, J. M. & Singh, G. Characterization of the mating-type locus (MAT) reveals a heterothallic mating system in Knightiella splachnirima. Lichenologist 49, 373–385 (2017).Article 

Google Scholar 
Czarnota, P. The lichen genus Micarea (Lecanorales, Ascomycota) in Poland. Pol. Bot. Stud. 23, 1–190 (2007).
Google Scholar 
Czarnota, P. & Guzow-Krzemińska, B. A phylogenetic study of the Micarea prasina group shows that Micarea micrococca includes three distinct lineages. Lichenologist 42, 7–21 (2010).Article 

Google Scholar 
Sérusiaux, E., Brand, A. M., Motiejūnaitè, J., Orange, A. & Coppins, B. J. Lecidea doliiformis belongs to Micarea, Catillaria alba to Biatora and Biatora lignimollis occurs in Western Europe. Bryologist 113, 333–344 (2010).Article 

Google Scholar 
van den Boom, P., Brand, A., Coppins, B. & Sérusiaux, E. Two new species in the Micarea prasina group from Western Europe. Lichenologist 49, 13–25 (2017).Article 

Google Scholar 
Guzow-Krzemińska, B., Czarnota, P., Łubek, A. & Kukwa, M. Micarea soralifera sp. nov., a new sorediate species in the M. prasina group. Lichenologist 48, 161–169 (2016).Article 

Google Scholar 
Guzow-Krzemińska, B. et al. Understanding the evolution of phenotypical characters in the Micarea prasina group (Pilocarpaceae) and descriptions of six new species within the group. MycoKeys 57, 1–30 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Kantvilas, G. & Coppins, B. J. Studies on Micarea in Australasia II. A synopsis of the genus in Tasmania, with the description of ten new species. Lichenologist 51, 431–481 (2019).Article 

Google Scholar 
Launis, A. & Myllys, L. Micarea fennica, a new lignicolous lichen species from Finland. Phytotaxa 409, 179–188 (2019).Article 

Google Scholar 
Launis, A., Pykälä, J., van den Boom, P., Sérusiaux, E. & Myllys, L. Four new epiphytic species in the Micarea prasina group from Europe. Lichenologist 51, 7–25 (2019).Article 

Google Scholar 
Launis, A. et al. Sharpening species boundaries in the Micarea prasina group, with a new circumscription of the type species M. prasina. Mycologia 111, 574–592 (2019).CAS 
PubMed 
Article 

Google Scholar 
van den Boom, P., Guzow-Krzemińska, B. & Kukwa, M. Two new Micarea species (Pilocarpaceae) from Western Europe. Plant Fungal Syst. 65, 189–199. https://doi.org/10.35535/pfsyst-2020-0014 (2020).Article 

Google Scholar 
Kantelinen, A., Hyvärinen, M., Kirika, P. & Myllys, L. Four new Micarea species from the montane cloud forests of Taita Hills, Kenya. Lichenologist 53, 81–94. https://doi.org/10.1017/S0024282920000511 (2021).Article 

Google Scholar 
Coppins, B. J. A taxonomic study of the lichen genus Micarea in Europe. Bull. Br. Mus. (Nat. Hist.) Bot. 11, 17–214 (1983).
Google Scholar 
Launis, A. & Myllys, L. Micarea byssacea new to North America and Micarea hedlundii new to Maine, Michigan and Quebec. Opusc. Philolichenum 13, 84–90 (2014).
Google Scholar 
Myllys, L. & Launis, A. Additions to the diversity of lichens and lichenicolous fungi living on decaying wood in Finland. Graphis Scr. 30, 78–87 (2018).
Google Scholar 
Yahr, R., Florence, A., Škaloud, P. & Voytsekhovich, A. Molecular and morphological diversity in photobionts associated with Micarea s. str. (Lecanorales, Ascomycota). Lichenologist 47, 403–414 (2015).Article 

Google Scholar 
Spribille, T., Thor, G., Bunnell, F. L., Goward, T. & Björk, C. R. Lichens on dead wood: Species-substrate relationships in the epiphytic lichen floras of the Pacific Northwest and Fennoscandia. Ecography 31, 741–750 (2008).Article 

Google Scholar 
Resl, P., Fernańdez-Mendoza, F., Mayrhofer, H. & Spribille, T. The evolution of fungal substrate specificity in a widespread group of crustose lichens. Proc. R. Soc. B 285, 20180640. https://doi.org/10.1098/rspb.2018.0640 (2018).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Stokland, J. N., Siitonen, J. & Jonsson, B. G. Biodiversity in Dead Wood 412 (Cambridge University Press, Cambridge, 2012).Book 

Google Scholar 
Russell, M. B., Woodall, C. W., Fraver, S. & D’Amato, A. W. Estimates of downed woody debris decay class transitions for forests across the eastern United States. Ecol. Model. 251, 22–31 (2013).Article 

Google Scholar 
Russell, M. B. et al. Residence times and decay rates of downed woody debris biomass/carbon in eastern US Forests. Ecosystems 17, 765–777 (2014).CAS 
Article 

Google Scholar 
Zoller, S., Lutzoni, F. & Scheidegger, C. Genetic variation within and among populations of the threatened lichen Lobaria pulmonaria in Switzerland and implications for its conservation. Mol. Ecol. 8, 2049–2059 (1999).CAS 
PubMed 
Article 

Google Scholar 
Honegger, R., Zippler, U., Gansner, H. & Scherrer, S. Mating systems in the genus Xanthoria (lichen forming Ascomycetes). Mycol. Res. 108, 480–488 (2004).CAS 
PubMed 
Article 

Google Scholar 
Honegger, R. & Zippler, U. Mating systems in representatives of the Parmeliaceae, Ramalinaceae and Physciaceae (Lecanoromycetes, lichen-forming ascomycetes). Mycol. Res. 11, 424–432 (2007).Article 
CAS 

Google Scholar 
Ament-Velásquez, S. L. et al. The plot thickens: Haploid and triploid-like thalli, hybridization, and biased mating Type Ratios in Letharia. Front. Fungal Biol. 2, 656386. https://doi.org/10.3389/ffunb.2021.656386 (2021).Article 

Google Scholar 
van den Boom, P. & Coppins, B. J. Micarea viridileprosa sp. nov., an overlooked lichen species from Western Europe. Lichenologist 33, 87–91 (2001).Article 

Google Scholar 
Simon, J.-C., Rispe, C. & Sunnucks, P. Ecology and evolution of sex in aphids. Trends Ecol. Evol. 17, 34–39 (2002).Article 

Google Scholar 
Silvertown, J. The evolutionary maintenance of sexual reproduction: Evidence from the ecological distribution of asexual reproduction in clonal plants. Int. J. Plant Sci. 169, 157–168 (2008).Article 

Google Scholar 
Gomez-Mestre, I., Pyron, R. A. & Wiens, J. J. Phylogenetic analyses reveal unexpected patterns in the evolution of reproductive modes in frogs. Evolution 66, 3687–3700. https://doi.org/10.1111/J.1558-5646.2012.01715.X (2012).Article 
PubMed 

Google Scholar 
Dańko, A., Schaible, R. & Dańko, M. J. Salinity effects on survival and reproduction of hydrozoan Eleutheria dichotoma. Estuaries Coasts 43, 360–374. https://doi.org/10.1007/s12237-019-00675-2 (2020).CAS 
Article 

Google Scholar 
Coppins, B. J. & Tønsberg, T. A new xanthone-containing Micarea from Northwest Europe and the Pacific Northwest of North America. Lichenologist 33, 93–96 (2001).Article 

Google Scholar 
Konoreva, L., Chesnokov, S., Kuznetsova, E. & Stepanchikova, I. Remarkable records of Micarea from the Russian Far East and significant extension of Micarea laeta and M. microareolata range. Botanica 25, 186–201 (2019).Article 

Google Scholar 
Weber, L., Printzen, C., Bässler, C. & Kantelinen, A. Seven Micarea (Pilocarpaceae) species new to Germany and notes on deficiently known species in the Bavarian forest. Herzogia 34, 5–17 (2021).Article 

Google Scholar 
van den Boom, P. Some interesting records of lichens and lichenicolous fungi from The Netherlands VI. Osten. Z. Pilzk. 12 (2003).Orange, A., James, P. W. & White, F. J. Microchemical Methods for the Identification of Lichens 101 (British Lichen Society, London, 2010).
Google Scholar 
Meyer, B. & Printzen, C. Proposal for a standardized nomenclature and characterization of insoluble lichen pigments. Lichenologist 32, 571–583 (2000).Article 

Google Scholar 
Culberson, C. F. & Kristinsson, H. D. A standardized method for the identification of lichen products. J. Chromatocraphy A 46, 85–93 (1970).CAS 
Article 

Google Scholar 
Myllys, L. et al. Phylogeny of the genus Bryoria. Lichenologist 43, 617–638 (2011).Article 

Google Scholar 
Myllys, L., Lohtander, K., Källersjö, M. & Tehler, A. Sequence insertion and ITS data provide congruent information in Roccella canariensis and R. tuberculata (Arthoniales, Euascomycetes) phylogeny. Mol. Phylogenetics Evol. 12, 295–309 (1999).CAS 
Article 

Google Scholar 
White, T. J., Bruns, T., Lee, S. & Taylor, J. W. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. In PCR Protocols: A Guide to the Methods and Applications (eds Innis, M. A. et al.) 315–322 (Academic Press, Cambridge, 1990).
Google Scholar 
Zoller, S., Scheidegger, C. & Sperisen, C. PCR primers for the amplification of mitochondrial small subunit ribosomal DNA of lichen-forming ascomycetes. Lichenologist 31, 511–516 (1999).Article 

Google Scholar 
Leavitt, S. D., Johnson, L. A., Goward, T. & Clair, L. L. S. Species delimitation in taxonomically difficult lichen-forming fungi: an example from morphologically and chemically diverse Xanthoparmelia (Parmeliaceae) in North America. Mol. Phylogenetics Evol. 60(3), 317–332 (2011).Article 

Google Scholar 
Schmitt, I. et al. New primers for promising single-copy genes in fungal phylogenetics and systematics. Persoonia 23, 35–40 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Edgar, R. C. MUSCLE: Multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797. https://doi.org/10.1093/nar/gkh340 (2004).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Kauff, F. & Lutzoni, F. Phylogeny of the Gyalectales and Ostropales (Ascomycota, Fungi): Among and within order relationships based on nuclear ribosomal RNA small and large subunits. Mol. Phylogenet. Evol. 25, 138–156 (2002).CAS 
PubMed 
Article 

Google Scholar 
Ronquist, F. & Huelsenbeck, J. P. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574 (2003).CAS 
PubMed 
Article 

Google Scholar 
Stamatakis, A. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Huelsenbeck, J. P., Larget, B. & Alfaro, M. E. Bayesian phylogenetic model selection using reversible jump markov chain monte carlo. Mol. Biol. Evol. 21, 1123–1133. https://doi.org/10.1093/molbev/msh123 (2004).CAS 
Article 
PubMed 

Google Scholar 
Rambaut, A., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior Summarization in bayesian phylogenetics using tracer 1.7. Syst. Biol. 67, 901–904 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Maddison, D. R. & Maddison, W. P. Mesquite: a modular system for evolutionary analysis. Version 3.40 http://mesquiteproject.org (2018).Pagel, M. Detecting correlated evolution on phylogenies: A general method for the comparative analysis of discrete characters. Proc. R. Soc. B. 255, 37–45 (1994).ADS 
Article 

Google Scholar  More

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.This is a summary of: Van der Stocken, T. et al. Mangrove dispersal disrupted by projected changes in global seawater density. Nat. Clim. Change https://doi.org/10.1038/s41558-022-01391-9 (2022). More

Sampling, DNA preparation and sequencingStockholmSamples LOW002, LOW003, LOW006, LOW007, LOW008 and PON012 were processed at the Archaeological Research Laboratory at Stockholm University, Sweden, following methods previously described8. In brief, this involved extracting DNA by incubating the bone powder for 24 h at 37 °C in 1.5 ml of digestion buffer (0.45 M EDTA (pH 8.0) and 0.25 mg ml–1 proteinase K), concentrating supernatant on Amicon Ultra-4 (30-kDa molecular weight cut-off (MWCO)) filter columns (MerckMillipore) and purifying on Qiagen MinElute columns. Double-stranded Illumina libraries were prepared using the protocol outlined in ref. 48, with the inclusion of USER enzyme and the modifications described in ref. 49.Samples 367, PDM100, Taimyr-1 and Yana-1 were processed at the Swedish Museum of Natural History in Stockholm, Sweden, following previously described methods8. In brief, this involved extracting DNA using a silica-based method with concentration on Vivaspin filters (Sartorius), according to a protocol optimized for recovery of ancient DNA50. Double-stranded Illumina libraries were prepared using the protocol outlined in ref. 48, with the inclusion of USER enzyme.Samples ALAS_024, VAL_033, ALAS_016, VAL_008, HMNH_007, HMNH_011, VAL_050, VAL_005, DS04, VAL_037, VAL_012, VAL_011, VAL_18A, IN18_016 and IN18_005 were processed at the Swedish Museum of Natural History in Stockholm, Sweden, following previously described methods for permafrost bone and tooth samples51. In brief, this involved DNA extraction using the methodology of ref. 52 and double-stranded Illumina library preparation as described in ref. 48, with dual unique indexes and the inclusion of USER enzyme. Between eight and ten separate PCR reactions with unique indexes were carried out for each sample to maximize library complexity. The libraries were sequenced alongside samples HOV4, AL2242, AL2370, AL2893, AL3272 and AL3284 across three Illumina NovaSeq 6000 lanes with an S4 100-bp paired-end set-up at SciLifeLab in Stockholm.PotsdamSamples JAL48, JAL65, JAL69, JAL358, AH574, AH575 and AH577 were processed at the University of Potsdam. Pre-amplification steps (DNA extraction and library preparation) were conducted in separated laboratory rooms specially equipped for the processing of ancient DNA. Amplification and post-amplification steps were performed in different laboratory rooms. DNA was extracted from bone powder (29–54 mg) following a protocol specially adapted to recover short DNA fragments52. Single-stranded double-indexed libraries were built from 20 µl of DNA extract according to the protocol in ref. 53. The libraries were sequenced on an HiSeq X platform at SciLifeLab in Stockholm.Tübingen/JenaSamples JK2174, JK2175, JK2179, JK2181, JK2183, TU144, TU148, TU839 and TU840 were processed at the University of Tübingen, with DNA extraction and pre-amplification steps undertaken in clean room facilities and post-amplification steps performed in a separate DNA laboratory. Both laboratories fulfil standards for work with ancient DNA54,55. All surfaces of tooth and bone samples were initially UV irradiated for 30 min, to minimize the potential risk of modern DNA contamination. Subsequently, DNA was extracted by applying a well-established guanidine silica-based protocol for ancient samples52. Illumina sequencing libraries were prepared by using 20 µl of DNA extract per library48; afterwards, dual barcodes (indexes) were chemically added to the prime ends of the libraries56. For the samples from Auneau (TU839 and TU840), five sequencing libraries each were prepared; for all other samples processed in Tübingen, three sequencing libraries each were prepared. To detect potential contamination of the chemicals, negative controls were conducted for extraction and library preparation. After preparation of the sequencing libraries, DNA concentration was measured with qPCR (Roche LightCycler) using corresponding primers48. The DNA concentration was given by the copy number of the DNA fragments in 1 µl of the sample.Amplification of the indexed sequencing libraries was performed using Herculase II Fusion under the following conditions: 1× Herculase II buffer, 0.4 µM IS5 primer and 0.4 µM IS6 primer48, Herculase II Fusion DNA polymerase (Agilent Technologies), 0.25 mM dNTPs (100 mM; 25 mM each dNTP) and 0.5–4 µl barcoded library as template in a total reaction volume of 100 µl. The applied amplification thermal profile was processed as follows: initial denaturation for 2 min at 95 °C; denaturation for 30 s at 95 °C, annealing for 30 s at 60 °C and elongation for 30 s at 72 °C for 3 to 20 cycles; and a final elongation step for 5 min at 72 °C. Thereafter, the amplified DNA was purified using a MinElute purification step and DNA was eluted in 20 µl TET. The concentration of the amplified DNA sequencing libraries was measured using a Bioanalyzer (Agilent Technologies) and a DNA1000 lab chip from Agilent Technologies.The sequencing libraries were sequenced on an Illumina HiSeq 4000 platform at the Max Planck Institute for Science of Human History in Jena. The samples from Auneau (TU839 and TU840) were paired-end sequenced applying 2 × 50 + 8 + 8 cycles. All other libraries prepared in Tübingen were single-end sequenced using 75 + 8 + 8 cycles.OxfordSamples AL2657, AL2541, AL2741, AL2744, AL3185, AL2350, CH1109, AL2370, AL3272 and AL3284 were processed at the dedicated ancient DNA facility at the PalaeoBARN laboratory at the University of Oxford, following methods described previously8. In brief, double-stranded libraries were constructed following the protocol in ref. 48. These libraries were sequenced on a HiSeq 2500 (AL2657, AL2541, AL2741, AL2744) or a HiSeq 4000 (AL3185, AL2350, CH1109) instrument at the Danish National Sequencing Center or on a NextSeq 550 instrument (AL2741) at the Natural History Museum of London. For samples AL2370, AL3272 and AL3284, between six and eight separate PCR reactions with unique indexes were carried out on their libraries and they were sequenced alongside samples HOV4, VAL_18A and IN18_016 on an Illumina NovaSeq 6000 lane with an S4 100-bp paired-end set-up at SciLifeLab in Stockholm.CopenhagenSamples CGG13, CGG17, CGG19, CGG20, CGG21, CGG25, CGG26, CGG27, CGG28, CGG34, Tumat1 and IRK were processed at the GLOBE Institute, University of Copenhagen. All pre-PCR work was performed in ancient DNA facilities following ancient DNA guidelines57. The details of extraction, library construction and sequencing for the samples with CGG codes are described in ref. 21, in relation to the publication of mitochondrial data from these specimens. The Tumat1 sample was processed following the exact same protocol. In brief, DNA extraction was performed using a buffer containing urea, EDTA and proteinase K50, double-stranded libraries were prepared with NEBNext DNA Sample Prep Master MixSet 2 (E6070S, New England Biolabs) and Illumina-specific adaptors48, and sequencing was performed on an Illumina HiSeq 2500 platform using 100-bp single-read chemistry. For the IRK sample, DNA was extracted from three subsamples and purified as described in ref. 21. The three DNA extracts and the purified pre-digest of one subsample were incorporated into double-stranded libraries following the BEST protocol58, with the modifications described in ref. 59, and sequenced on a BGISEQ-500 platform using 100-bp single-read chemistry.Santa CruzSamples SC19.MCJ017, SC19.MCJ015, SC19.MCJ010 and SC19.MCJ014 were processed at the UCSC Paleogenomics Lab and were provided by the Yukon Government Paleontology program. All pre-PCR work was performed in a dedicated ancient DNA facility at the University of California, Santa Cruz, following standard ancient DNA methods60. Subsamples (250–350 mg) were sent to the UCI KECK AMS facility for radiocarbon dating, and the remaining amounts were powdered in a Retsch MM400 for extraction. For each sample, ~100 mg of powder was treated with a 0.5% sodium hypochlorite solution before extraction to remove surface contaminants61 and then combined with 1 ml lysis buffer for extraction, following the protocol in ref. 52. Samples were processed in parallel with a negative control. We quantified the extracts using a Qubit 1× dsDNA HS Assay kit (Q33231) before preparing libraries. We prepared single-stranded libraries following the protocol in ref. 62 and amplified the libraries for 9–16 cycles as informed by qPCR. After amplification, we cleaned the libraries using a 1.2× SPRI bead solution and pooled them to an equimolar ratio for in-house shallow quality-control sequencing on a NextSeq 550 paired-end 75-bp run. We then sent the libraries to Fulgent Genetics for deeper sequencing on two paired-end 150-bp lanes on a HiSeq X instrument.ViennaSample HOV4 was processed at the Department of Anthropology, University of Vienna. The sample is a canine tooth, which after sequencing was determined to derive from a dhole (Cuon alpinus). DNA was extracted from its cementum using the methods described in ref. 63 with a modified incubation time of ~18 h. The library was prepared according to the protocol in ref. 48 with the modifications from ref. 64. Five separate PCR reactions with unique indexes were carried out on the library and were sequenced alongside samples VAL_18A, IN18_016, AL2242, AL2370, AL2893, AL3272 and AL3284 on an Illumina NovaSeq 6000 lane with an S4 100-bp paired-end set-up at SciLifeLab in Stockholm.An overview of all samples and their associated metadata is available in Supplementary Data 1.Genome sequence data processingFor paired-end data, read pairs were merged and adaptors were trimmed using SeqPrep (https://github.com/jstjohn/SeqPrep), discarding reads that could not be successfully merged. Reads were mapped to the dog reference genome canFam3.1 using BWA aln (v.0.7.17)65 with permissive parameters, including a disabled seed (-l 16500 -n 0.01 -o 2). Duplicates were removed by keeping only one read from any set of reads that had the same orientation, length and start and end coordinates. For sample Taimyr-1, previously published data13 were merged with newly generated data. Data from samples processed in Copenhagen were processed as described previously66 except that they were also mapped to canFam3.1. Post-mortem damage was quantified using PMDtools (v0.60)67 with the ‘–first’ and ‘–CpG’ arguments.Genotyping and integration with previously published genomesTo construct a comparative dataset for population genetic analyses, we started from a published variant call set compiling 722 modern dog, wolf and other canid genomes from multiple previous studies (NCBI BioProject accession PRJNA448733)40. To this, we added additional modern whole genomes from other studies: 4 African golden wolves and 15 Nigerian village dogs (Genome Sequence Archive (http://gsa.big.ac.cn/), accession PRJCA000335)68, 12 Scandinavian wolves (European Nucleotide Archive accession PRJEB20635)69, 9 North American wolves and coyotes (European Nucleotide Archive accession PRJNA496590)25 and 8 other canids (African hunting dog, dhole, Ethiopian wolf, golden jackal, Middle Eastern grey wolves) (European Nucleotide Archive accession PRJNA494815)22. Reads from these genomes were mapped to the dog reference genome using bwa mem (version 0.7.15)70, marked for duplicates using Picard Tools (v2.21.4) (http://broadinstitute.github.io/picard), genotyped at the sites present in the above dataset using GATK HaplotypeCaller (v3.6)71 with the ‘-gt_mode GENOTYPE_GIVEN_ALLELES’ argument and then merged into the dataset using bcftools merge (http://www.htslib.org/). The following filters were then applied to sites and genotypes across the full dataset: sites with excess heterozygosity (bcftools fill-tags ‘ExcHet’ P value 5.8×. For divergence time analyses, haploid X chromosomes from two different male genomes were combined and the point at which the inferred effective population size for this ‘pseudodiploid’ chromosome increased sharply upwards was taken to correspond to a population divergence. Results were scaled using a mutation rate of 0.4 × 10−8 mutations per site per generation13,87 (with a 25% lower rate for X-chromosome analyses) and a mean generational interval of 3 years13. For effective population size inferences, transition variants were ignored and results were scaled using a transversions-only mutation rate inferred from results on modern genomes. For more details on the MSMC2 analyses, see Supplementary Information section 3.Selection analysesSelection analysis was performed using PLINK (v1.90b5.2)88. This analysis used the 72 ancient wolf genomes and 68 modern wolf genomes (with the latter including a historical Japanese wolf genome73 treated as ancient for analysis purposes, with its age set to 200 bp). A list of the genomes used for this analysis is available in Supplementary Data 2 (“Used for selection scan” column). All SNPs, not only transversions, were used for this analysis. The age of each wolf was set as the phenotype, with values of 0 for modern wolves, and the ‘–linear’ argument was used to test for an association between SNP genotypes and age, also applying the ‘–adjust’ argument to correct P values using genomic control. The application of genomic control34 here aimed to use the magnitude of temporal allele frequency variance observed across the genome to account for what was observed from genetic drift alone given wolf demographic history. Only results for the following sets of sites were retained and included in the Manhattan plot: sites where at least 40 ancient genomes had a genotype call, sites with a minor allele frequency among the ancient wolves of ≥5% and sites that had at least 7 neighbouring sites within a 50-kb window with a P value that was at least 90% as large (on a log10 scale) as the P value of the site itself. The last ‘neighbourhood filter’ aimed to reduce false positives by requiring similar evidence across multiple nearby sites. As a P-value significance cut-off to correct for the genome-wide testing, we used 5 × 10−8, which is commonly used in genome-wide association studies in humans and also in dogs89. We excluded 15 regions where only a single variant reached significance. A detailed table with the 24 detected regions is available in Supplementary Data 3. To test the robustness of this analysis to false positives arising from genetic drift alone, we applied the same analysis to data from neutral coalescent simulations generated using ms90 and found no false positives. For more details, see Supplementary Information section 4.Ancestry modelling with qpAdm and qpWaveWe used the qpAdm and qpWave methods43 from ADMIXTOOLS (v5.0)84 to test ancestry models for wolf and dog targets postdating 23 ka. For the primary analyses, we used the following set of candidate source populations (age estimate in brackets, years bp): Armenia_Hovk1.HOV4 (ancient dhole), Siberia_UlakhanSular.LOW008 (70,772), Germany_Aufhausener.AH575 (57,233), Siberia_BungeToll.CGG29 (48,210), Germany_HohleFels.JK2183 (32,366), Siberia_BelayaGora.IN18_016 (32,020), Yukon_QuartzCreek.SC19.MCJ010 (29,943), Altai_Razboinichya.AL2744 (28,345), Siberia_BelayaGora.IN18_005 (18,148) and Germany_HohleFels.JK2179 (13,229). We used a rotating approach in which, for each target, we tested all possible one-, two- and three-source models that could be enumerated from the above set. Individuals from the set that were not used as a source in a given model served as thereference set (or the ‘right’ population in the qpAdm framework). This means that, in every model, each of the above individuals was always either in the source list or in the reference list. We ranked models on the basis of their P values, but prioritized models with fewer sources using a P-value threshold of 0.01: if a simpler model (meaning a model with fewer sources) had a P value above this threshold, it ranked above a more complex model (meaning a model with more sources) regardless of the P value of the latter. We also failed models with inferred ancestry proportions larger than 1.1 or smaller than −0.1. For single-source models, qpWave was run instead of qpAdm. Both programs were run with the ‘allsnps: YES’ option (without this option, there was very little power to reject models). We describe ancestry assigned to the ancient dhole source (Armenia_Hovk1.HOV4) as ‘unsampled’ ancestry; note that this does not imply that such ancestry is of non-wolf origin, only that it is not represented by (that is, diverged early from and lacks shared genetic drift with) the ancient wolf genomes in the reference set.To test whether any post-23 ka or modern wolf genome available might be a good proxy for the western Eurasian wolf-related ancestry identified in Near Eastern and African dogs, we added the 9,500-year-old Zhokhov dog17 to the rotating set of candidate source populations. Chosen for its high coverage, early date and easterly location, this makes the assumption that the Zhokhov dog is a good representative for the eastern dog ancestry component. Using the African Basenji dog as a target, models involving the Zhokhov dog plus another given wolf thus allowed us to test whether that wolf was a good match for the additional component of ancestry. For more details on the qpAdm and qpWave analyses, see Supplementary Information sections 2 (wolf targets) and  5 (dog targets).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this paper. More