Ecology

Physico-chemical conditionsSampling was performed at 6 stations representing the physical and chemical characteristics of the investigated area (Supplementary Table S1). Thermohaline properties were the result of horizontal advection of above-average salinities driven by a North Ionian cyclonic gyre controlled by the Adriatic Ionian Bimodal Oscillating System46. September and the whole summer of 2016 was characterized by extremely high temperatures, and precipitation in the climatologic expected range. A cyclone with a cold front followed by a strong Bora wind passed over the Adriatic a week before the cruise, in the period between the 16th and 20th of September 2016. Heat and mass exchange in the air-sea boundary layer were responsible for the characteristic vertical thermohaline profiles measured in late summer. Over the investigated area, the mixed layer depth located between 20 and 25 m was horizontally homogenous. The coldest water mass (temperature 12.94 °C, salinity 38.68) was located at the bottom of Jabuka Pit.Abundance of bacteria, autotrophic picoplankton and AAPBacterial abundances ranged between 0.05 and 0.46 × 106 cell mL−1 in all three areas, with a slightly higher average value in Jabuka Pit (0.31 × 106 cell mL−1). The bacterial abundances were the highest in the upper layers down to the 50 m deep layer and showed a decreasing trend towards the bottom (Supplementary Table S2). The portion of HNA bacteria ranged from 37.8 to 73.12% (on average 51.27%), with the prevalence of HNA over the LNA group below the epipelagic layer.Marine Synechococcus dominated the autotrophic picoplankton community with abundances ranging from 0.08 to 23.86 × 103 cell mL−1. The presence of Prochlorococcus cells was also detected in all samples in a range from a few cells to 1.33 × 103 cell mL−1. Picoeukaryotes also showed a similar range from a few cells to 0.83 × 103 cell mL−1. The highest abundances of picophytoplankton were measured in the upper 50 m, with the exception of the Palagruža Sill (PS) area, where an increase in abundance was observed at 100 m depth. Bacterial production ranged from 0.2 × 104 to 0.36 × 104 cell mL−1 h−1, with increased values in the shallow layers and a mostly uniform vertical distribution in the water column (Supplementary Table S2).AAP bacteria abundance ranged from 0.9 × 103 to 22.3 × 103 cell mL−1, thus constituting 0.42% to 6.83% of the bacteria. Their highest average contribution was observed in the South Adriatic Pit (4.11%), while on the vertical scale, their highest contribution was observed in the upper 20 m of the seawater (see Supplementary Table S2).Relationship between the picoplankton community and environmental parametersBased on biological characteristics (total prokaryotes, Synechococcus, Prochlorococcus, picoeukaryotes, heterotrophic nanoflagellates, aerobic anoxygenic phototrophs, high and low nucleic acid bacteria, bacterial production), we distinguished five picoplanktonic clusters (PIC-BMUs) and then searched for explanations of the observed patterns (Fig. 2A,B). The mean values of biological and physico-chemical parameters for each cluster are shown in Table 1.Figure 2(A) Bar plot representation of biotic (black) and abiotic (grey) parameters for neural gas best-matching units (picoplankton-PIC-BMUs) with relative frequency appearance for each neuron. TP-total prokaryotes, SYN-Synechococcus, PROCHL-Prochlorococcus, PE-picoeukaryotes, HNF-heterotrophic nanoflagellates, AAP-aerobic anoxygenic phototrophs, AAP%-portion of AAP, HNA% percentage of high nucleic acid content bacteria, LNA%-percentage of low nucleic acid content bacteria-LNA%, BP-bacterial production. (B) Water column distribution of Neural gas best-matching units (BMU, labels with numbers, and stained with a different colour for clearance, coloured non-labelled squares shows clarity) for measuring stations (SAP1-3, PS1-2 and JP1). The software MATLAB. version 7.10.0 (R2018). Natick, Massachusetts: The MathWorks Inc. (2018) (https://www.mathworks.com/) was used to generate the figure.Full size imageTable 1 Characteristics of biological (abundances of total prokaryotes-TP, Synechococcus-SYN, Prochlorococcus-PROCHL, picoeukaryotes-PE, heterotrophic nanoflagellates-HNF, aerobic anoxygenic phototrophs(AAP); contributions (%) of AAP, High nucleic acid content bacteria-HNA and Low nucleic acid content bacteria-LNA%; and bacterial production-BP) and environmental factors in the sampling terms assigned to the neural gas clusters.Full size tablePIC-BMU1 described a very rare pattern, found in only two samples from 10 m depth in Palaguža Sill and Jabuka Pit. They were characterised by the highest abundances of total prokaryotes with a dominance of HNA and elevated AAP abundance. These samples were unique in terms of hydrological parameters, as they represented an N-limited environment (TIN  More

Aivelo T, Norberg A (2018) Parasite-microbiota interactions potentially affect intestinal communities in wild mammals. J Anim Ecol 87:438–447PubMed 
Article 
PubMed Central 

Google Scholar 
Alberdi A, Aizpurua O, Bohmann K, Zepeda-Mendoza ML, Gilbert MTP (2016) Do vertebrate gut metagenomes confer rapid ecological adaptation? Trends Ecol Evol 31:689–699PubMed 
Article 
PubMed Central 

Google Scholar 
Amaral WZ, Lubach GR, Proctor A, Lyte M, Phillips GJ, Coe CL (2017) Social influences on Prevotella and the gut microbiome of young monkeys. Psychosom Med 79:888–897PubMed 
PubMed Central 
Article 

Google Scholar 
Amato KR, Sanders GJ, Song SJ, Nute M, Metcalf JL, Thompson LR et al. (2019) Evolutionary trends in host physiology outweigh dietary niche in structuring primate gut microbiomes. ISME J 13:576–587CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol 57:289–300
Google Scholar 
Björk JR, Dasari M, Grieneisen L, Archie EA (2019) Primate microbiomes over time: longitudinal answers to standing questions in microbiome research. Am J Primatol 81:e22970PubMed 
PubMed Central 
Article 

Google Scholar 
Brooks JW (2016) Postmortem changes in animal carcasses and estimation of the postmortem interval. Vet Pathol 53:929–940CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP (2016) DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods 13:581–583CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Callahan BJ, Wong J, Heiner C, Oh S, Theriot CM, Gulati AS et al. (2019) High-throughput amplicon sequencing of the full-length 16S rRNA gene with single-nucleotide resolution. Nucleic Acids Res 47:e103CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Clayton JB, Vangay P, Huang H, Ward T, Hillmann BM, Al-Ghalith GA et al. (2016) Captivity humanizes the primate microbiome. Proc Natl Acad Sci U S A 113:10376–10381CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Cryan JF, Dinan TG (2012) Mind-altering microorganisms: the impact of the gut microbiota on brain and behaviour. Nat Rev Neurosci 13:701–712CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
De Filippo C, Cavalieri D, Di Paola M, Ramazzotti M, Poullet JB, Massart S et al. (2010) Impact of diet in shaping gut microbiota revealed by a comparative study in children from Europe and rural Africa. Proc Natl Acad Sci U S A 107:14691–14696PubMed 
PubMed Central 
Article 

Google Scholar 
Dill-McFarland KA, Neil KL, Zeng A, Sprenger RJ, Kurtz CC, Suen G et al. (2014) Hibernation alters the diversity and composition of mucosa-associated bacteria while enhancing antimicrobial defence in the gut of 13-lined ground squirrels. Mol Ecol 23:4658–4669CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Donaldson GP, Lee SM, Mazmanian SK (2016) Gut biogeography of the bacterial microbiota. Nat Rev Microbiol 14:20–32CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Dubois S, Fenwick N, Ryan EA, Baker L, Baker SE, Beausoleil NJ et al. (2017) International consensus principles for ethical wildlife control. Conserv Biol J Soc Conserv Biol 31:753–760Article 

Google Scholar 
Earl JP, Adappa ND, Krol J, Bhat AS, Balashov S, Ehrlich RL et al. (2018) Species-level bacterial community profiling of the healthy sinonasal microbiome using Pacific Biosciences sequencing of full-length 16S rRNA genes. Microbiome 6:190PubMed 
PubMed Central 
Article 

Google Scholar 
Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R (2011) UCHIME improves sensitivity and speed of chimera detection. Bioinformatics 27:2194–2200CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ericsson AC, Johnson PJ, Lopes MA, Perry SC, Lanter HR (2016) A microbiological map of the healthy equine gastrointestinal tract. PLoS ONE 11:e0166523PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
García-Amado MA, Michelangeli F, Gueneau P, Perez ME, Domínguez-Bello MG (2007) Bacterial detoxification of saponins in the crop of the avian foregut fermenter Opisthocomus hoazin. J Anim Feed Sci 16:82–85Article 

Google Scholar 
Gomez A, Petrzelkova KJ, Burns MB, Yeoman CJ, Amato KR, Vlckova K et al. (2016) Gut microbiome of coexisting BaAka pygmies and Bantu reflects gradients of traditional subsistence patterns. Cell Rep 14:2142–2153CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Gomez A, Petrzelkova K, Yeoman CJ, Vlckova K, Mrázek J, Koppova I et al. (2015) Gut microbiome composition and metabolomic profiles of wild western lowland gorillas (Gorilla gorilla gorilla) reflect host ecology. Mol Ecol 24:2551–2565CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Gorvitovskaia A, Holmes SP, Huse SM (2016) Interpreting Prevotella and bacteroides as biomarkers of diet and lifestyle. Microbiome 4:15PubMed 
PubMed Central 
Article 

Google Scholar 
Gregorová S, Forejt J (2000) PWD/Ph and PWK/Ph inbred mouse strains of Mus m. musculus subspecies–a valuable resource of phenotypic variations and genomic polymorphisms. Folia Biol 46:31–41
Google Scholar 
Gu S, Chen D, Zhang J-N, Lv X, Wang K, Duan L-P et al. (2013) Bacterial community mapping of the mouse gastrointestinal tract. PLoS ONE 8:e74957Heimesaat MM, Boelke S, Fischer A, Haag L-M, Loddenkemper C, Kühl AA et al. (2012) Comprehensive postmortem analyses of intestinal microbiota changes and bacterial translocation in human flora associated mice. PloS ONE 7:e40758CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hird SM (2017) Evolutionary biology needs wild microbiomes. Front Microbiol 8:725Iljazovic A, Roy U, Gálvez EJC, Lesker TR, Zhao B, Gronow A et al. (2020) Perturbation of the gut microbiome by Prevotella spp. enhances host susceptibility to mucosal inflammation. Mucosal Immunol 14:113–124PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Ingala MR, Simmons NB, Wultsch C, Krampis K, Speer KA, Perkins SL (2018) Comparing microbiome sampling methods in a wild mammal: fecal and intestinal samples record different signals of host ecology, evolution. Front Microbiol 9:803Karasov WH, Douglas AE (2013) Comparative digestive physiology. Compr Physiol 3:741–783PubMed 
PubMed Central 

Google Scholar 
Kartzinel TR, Hsing JC, Musili PM, Brown BRP, Pringle RM (2019) Covariation of diet and gut microbiome in African megafauna. Proc Natl Acad Sci U S A 116:23588–23593CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kohl KD, Dearing MD (2016) The woodrat gut microbiota as an experimental system for understanding microbial metabolism of dietary toxins. Front Microbiol 7:1165Kohl KD, Luong K, Dearing MD (2015) Validating the use of trap-collected feces for studying the gut microbiota of a small mammal (Neotoma lepida). J Mammal 96:90–93Article 

Google Scholar 
Kohl KD, Varner J, Wilkening JL, Dearing MD (2018) Gut microbial communities of American pikas (Ochotona princeps): Evidence for phylosymbiosis and adaptations to novel diets. J Anim Ecol 87:323–330PubMed 
Article 
PubMed Central 

Google Scholar 
Kreisinger J, Bastien G, Hauffe HC, Marchesi J, Perkins SE (2015) Interactions between multiple helminths and the gut microbiota in wild rodents. Philos Trans R Soc B Biol Sci 370:20140295Kreisinger J, Čížková D, Vohánka J, Piálek J (2014) Gastrointestinal microbiota of wild and inbred individuals of two house mouse subspecies assessed using high-throughput parallel pyrosequencing. Mol Ecol 23:5048–5060CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Kreisinger J, Kropáčková L, Petrželková A, Adámková M, Tomášek O, Martin J-F et al. (2017) Temporal stability and the effect of transgenerational transfer on fecal microbiota structure in a long distance migratory bird. Front Microbiol 8:50Laukaitis CM, Critser ES, Karn RC (1997) Salivary androgen-binding protein (ABP) mediates sexual isolation in Mus musculus. Evol Int J Org Evol 51:2000–2005CAS 
Article 

Google Scholar 
Lawrence K, Lam K, Morgun A, Shulzhenko NLöhr C (2019) Effect of temperature and time on the thanatomicrobiome of the cecum, ileum, kidney, and lung of domestic rabbits. J Vet Diagn Invest 31. https://doi.org/10.1177/1040638719828412Legendre P, Anderson MJ (1999) Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments. Ecol Monogr 69:1–24Article 

Google Scholar 
Li D, Chen H, Mao B, Yang Q, Zhao J, Gu Z et al. (2017) Microbial biogeography and core microbiota of the rat digestive tract. Sci Rep 7:45840PubMed Central 
Article 
CAS 

Google Scholar 
Maslanik T, Tannura K, Mahaffey L, Loughridge AB, Benninson L, Ursell L et al. (2012) Commensal bacteria and MAMPs are necessary for stress-induced increases in IL-1β and IL-18 but not IL-6, IL-10 or MCP-1. PLoS ONE 7:e50636CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Matsuo Y, Komiya S, Yasumizu Y, Yasuoka Y, Mizushima K, Takagi T et al. (2020) Full-length 16S rRNA gene amplicon analysis of human gut microbiota using MinIONTM nanopore sequencing confers species-level resolution. bioRxiv. https://doi.org/10.1101/2020.05.06.078147McKenzie VJ, Song SJ, Delsuc F, Prest TL, Oliverio AM, Korpita TM et al. (2017) The effects of captivity on the mammalian gut microbiome. Integr Comp Biol 57:690–704PubMed 
PubMed Central 
Article 

Google Scholar 
McMurdie PJ, Holmes S (2014) Waste not, want not: why rarefying microbiome data is inadmissible. PLOS Comput Biol 10:e1003531PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Menke S, Meier M, Sommer S (2015) Shifts in the gut microbiome observed in wildlife faecal samples exposed to natural weather conditions: lessons from time-series analyses using next-generation sequencing for application in field studies. Methods Ecol Evol 6:1080–1087Article 

Google Scholar 
Miller AW, Oakeson KF, Dale C, Dearing MD (2016) Microbial community transplant results in increased and long-term oxalate degradation. Micro Ecol 72:470–478CAS 
Article 

Google Scholar 
Pafčo B, Čížková D, Kreisinger J, Hasegawa H, Vallo P, Shutt K et al. (2018) Metabarcoding analysis of strongylid nematode diversity in two sympatric primate species. Sci Rep 8:5933Palm NW, de Zoete MR, Cullen TW, Barry NA, Stefanowski J, Hao L et al. (2014) Immunoglobulin A coating identifies colitogenic bacteria in inflammatory bowel disease. Cell 158:1000–1010CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Pechal JL, Schmidt CJ, Jordan HR, Benbow ME (2018) A large-scale survey of the postmortem human microbiome, and its potential to provide insight into the living health condition. Sci Rep 8:5724PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Pollock J, Glendinning L, Wisedchanwet T, Watson M (2018) The madness of microbiome: attempting to find consensus “best practice” for 16S microbiome studies. Appl Environ Microbiol 84:e02627–17PubMed 
PubMed Central 
Article 

Google Scholar 
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P et al. (2013) The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res 41:D590–D596CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
R Core Team (2018) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/
Google Scholar 
Rosshart SP, Vassallo BG, Angeletti D, Hutchinson DS, Morgan AP, Takeda K et al. (2017) Wild mouse gut microbiota promotes host fitness and improves disease resistance. Cell 171:1015–1028.e13CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Round JL, Mazmanian SK (2009) The gut microbiome shapes intestinal immune responses during health and disease. Nat Rev Immunol 9:313–323CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Scher JU, Sczesnak A, Longman RS, Segata N, Ubeda C, Bielski C et al. (2013) Expansion of intestinal Prevotella copri correlates with enhanced susceptibility to arthritis. eLife 2:e01202PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Sommer F, Ståhlman M, Ilkayeva O, Arnemo JM, Kindberg J, Josefsson J et al. (2016) The gut microbiota modulates energy metabolism in the hibernating brown bear Ursus arctos. Cell Rep 14:1655–1661CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Stalder GL, Pinior B, Zwirzitz B, Loncaric I, Jakupović D, Vetter SG et al. (2019) Gut microbiota of the European Brown Hare (Lepus europaeus). Sci Rep 9:2738CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Stanley D, Geier MS, Chen H, Hughes RJ, Moore RJ (2015) Comparison of fecal and cecal microbiotas reveals qualitative similarities but quantitative differences. BMC Microbiol 15:51PubMed 
PubMed Central 
Article 

Google Scholar 
Stearns JC, Lynch MDJ, Senadheera DB, Tenenbaum HC, Goldberg MB, Cvitkovitch DG et al. (2011) Bacterial biogeography of the human digestive tract. Sci Rep 1:170CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Stothart MR, Palme R, Newman AEM (2019) It’s what’s on the inside that counts: stress physiology and the bacterial microbiome of a wild urban mammal. Proc R Soc B Biol Sci 286:20192111Article 

Google Scholar 
Suzuki TA, Martins FM, Nachman MW (2019) Altitudinal variation of the gut microbiota in wild house mice. Mol Ecol 28:2378–2390CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Suzuki TA, Nachman MW (2016) Spatial heterogeneity of gut microbial composition along the gastrointestinal tract in natural populations of house mice (EG Zoetendal, Ed.). PLoS ONE 11:e0163720PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Tanca A, Manghina V, Fraumene C, Palomba A, Abbondio M, Deligios M et al. (2017) Metaproteogenomics reveals taxonomic and functional changes between cecal and fecal microbiota in mouse. Front Microbiol 8:391Tang Q, Jin G, Wang G, Liu T, Liu X, Wang B et al. (2020) Current sampling methods for gut microbiota: a call for more precise devices. Front Cell Infect Microbiol 10:151Tang W, Zhu G, Shi Q, Yang S, Ma T, Mishra SK et al. (2019) Characterizing the microbiota in gastrointestinal tract segments of Rhabdophis subminiatus: dynamic changes and functional predictions. MicrobiologyOpen 8:e789Trevelline BK, Fontaine SS, Hartup BK, Kohl KD (2019) Conservation biology needs a microbial renaissance: a call for the consideration of host-associated microbiota in wildlife management practices. Proc R Soc B Biol Sci 286:20182448Article 

Google Scholar 
Tuomisto S, Karhunen PJ, Pessi T (2013) Time-dependent post mortem changes in the composition of intestinal bacteria using real-time quantitative PCR. Gut Pathog 5:35Turnbaugh PJ, Ley RE, Mahowald MA, Magrini V, Mardis ER, Gordon JI (2006) An obesity-associated gut microbiome with increased capacity for energy harvest. Nature 444:1027–1031PubMed 
Article 
PubMed Central 

Google Scholar 
Vasemägi A, Visse M, Kisand V (2017) Effect of Environmental Factors and an Emerging Parasitic Disease on Gut Microbiome of Wild Salmonid Fish. mSphere 2:e00418–17PubMed 
PubMed Central 
Article 

Google Scholar 
Videvall E, Strandh M, Engelbrecht A, Cloete S, Cornwallis C (2017) Measuring the gut microbiome in birds: Comparison of faecal and cloacal sampling. Mol Ecol Resour 18:424–434PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
Vlčková K, Shutt-Phillips K, Heistermann M, Pafčo B, Petrželková KJ, Todd A et al. (2018) Impact of stress on the gut microbiome of free-ranging western lowland gorillas. Microbiol Read Engl 164:40–44Article 
CAS 

Google Scholar 
Wang Q, Garrity GM, Tiedje JM, Cole JR (2007) Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microbiol 73:5261–5267CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Wang J, Linnenbrink M, Künzel S, Fernandes R, Nadeau M-J, Rosenstiel P et al. (2014) Dietary history contributes to enterotype-like clustering and functional metagenomic content in the intestinal microbiome of wild mice. Proc Natl Acad Sci U S A 111:E2703–2710CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Warne RW, Kirschman L, Zeglin L (2017) Manipulation of gut microbiota reveals shifting community structure shaped by host developmental windows in amphibian larvae. Integr Comp Biol 57:786–794PubMed 
Article 
PubMed Central 

Google Scholar 
Weldon L, Abolins S, Lenzi L, Bourne C, Riley EM, Viney M (2015) The gut microbiota of wild mice. PLoS ONE 10:e0134643PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Wu GD, Chen J, Hoffmann C, Bittinger K, Chen Y-Y, Keilbaugh SA et al. (2011) Linking long-term dietary patterns with gut microbial enterotypes. Science 334:105–108CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Yan W, Sun C, Zheng J, Wen C, Ji C, Zhang D et al. (2019) Efficacy of fecal sampling as a gut proxy in the study of chicken gut microbiota. Front Microbiol 10:2126Yasuda K, Oh K, Ren B, Tickle TL, Franzosa EA, Wachtman LM et al. (2015) Biogeography of the intestinal mucosal and lumenal microbiome in the rhesus macaque. Cell Host Microbe 17:385–391CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Zemanova MA (2019) Poor implementation of non-invasive sampling in wildlife genetics studies. Rethink Ecol 4:119–132Article 

Google Scholar 
Zemanova MA (2020) Towards more compassionate wildlife research through the 3Rs principles: moving from invasive to non-invasive methods. Wildl Biol 2020. https://doi.org/10.2981/wlb.00607Zhao W, Wang Y, Liu S, Huang J, Zhai Z, He C et al. (2015) The dynamic distribution of porcine microbiota across different ages and gastrointestinal tract segments. PLoS ONE 10:e0117441PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Zilber-Rosenberg I, Rosenberg E (2008) Role of microorganisms in the evolution of animals and plants: the hologenome theory of evolution. FEMS Microbiol Rev 32:723–735CAS 
PubMed 
Article 
PubMed Central 

Google Scholar  More

1.Kimura, T. Systemic alopecia resulting from hyperadrenocorticism in a Japanese monkey. Lab. Primate Newsl. 47, 5–9 (2008).
Google Scholar 
2.Novak, M. A. et al. Assessing significant ( > 30%) alopecia as a possible biomarker for stress in captive rhesus monkeys (Macaca mulatta). Am. J. Primatol. 79, e22547 (2017).Article 
CAS 

Google Scholar 
3.Lutz, C. K., Menard, M. T., Rosenberg, K., Meyer, J. S. & Novak, M. A. Alopecia in rhesus macaques (Macaca mulatta): Association with pregnancy and chronic stress. J. Med. Primatol. 48, 251–256. https://doi.org/10.1111/jmp.12419 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
4.Steinmetz, H. W., Kaumanns, W., Dix, I., Neimeier, K.-A. & Kaup, F.-J. Dermatologic investigation of alopecia in rhesus macaques (Macaca mulatta). J. Zoo. Wildl. Med. 36, 229–239 (2005).PubMed 
Article 

Google Scholar 
5.Lynch, M., Kirkwood, R., Mitchell, A., Duignan, P. & Arnould, J. P. Y. Prevalence and significance of an alopecia syndrome in Australian fur seals (Arctocephalus pusillus doriferus). J. Mammal. 92, 342–351 (2011).Article 

Google Scholar 
6.Atwood, T. et al. Prevalence and spatio-temporal variation of an alopecia syndrome in polar bears (Ursus maritimus) of the southern Beaufort Sea. J. Wildl. Dis. 51, 48–59 (2015).PubMed 
Article 

Google Scholar 
7.McCoy, R. H., Murphie, S. L., Szykman Gunther, M. & Murphie, B. L. Influence of hair loss syndrome on black-tailed deer fawn survival. J. Wildl. Manag. 78, 1177–1188. https://doi.org/10.1002/jwmg.772 (2014).Article 

Google Scholar 
8.Novak, M. A. et al. Hair loss and hypothalamic–pituitary–adrenocortical axis activity in captive rhesus macaques (Macaca mulatta). J. Am. Assoc. Lab. Anim. Sci. 53, 261–266 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
9.Lutz, C. K., Coleman, K., Worlein, J. & Novak, M. A. Hair loss and hair-pulling in rhesus macaques (Macaca mulatta). J. Am. Assoc. Lab. Anim. Sci. 52, 454–457 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
10.Swenerton, H. & Hurley, L. S. Zinc deficiency in rhesus and bonnet monkeys, including effects on reproduction. J. Nutr. 110, 575–583 (1980).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Lair, S., Crawshaw, G. J., Mehren, K. G. & Perrone, M. A. Diagnosis of hypothyroidism in a western lowland gorilla (Gorilla gorilla gorilla) using human thyroid-stimulating hormone assay. J. Zoo. Wildl. Med. 30, 537–540 (1999).CAS 
PubMed 
PubMed Central 

Google Scholar 
12.Beardi, B. et al. Alopecia areata in a rhesus monkey (Macaca mulatta). J. Med. Primatol. 36, 124–130. https://doi.org/10.1111/j.1600-0684.2007.00212.x (2007).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
13.Ovadia, S., Wilson, S. R. & Zeiss, C. J. Successful cyclosporine treatment for atopic dermatitis in a rhesus macaque (Macaca mulatta). Comp. Med. 55, 192–196 (2005).CAS 
PubMed 
PubMed Central 

Google Scholar 
14.Novak, M. A. & Meyer, J. S. Alopecia: Possible causes and treatments, particularly in captive nonhuman primates. Comp. Med. 59, 18–26 (2009).CAS 
PubMed 
PubMed Central 

Google Scholar 
15.Hadshiew, I. M., Foitzik, K., Arck, P. C. & Paus, R. Burden of hair loss: Stress and the underestimated psychosocial impact of telogen effluvium and androgenetic alopecia. J. Investig. Dermatol. 123, 455–457 (2004).CAS 
PubMed 
Article 

Google Scholar 
16.Charmandari, E., Tsigos, C. & Chrousos, G. Endocrinology of the stress response. Annu. Rev. Physiol. 67, 259–284 (2005).CAS 
PubMed 
Article 

Google Scholar 
17.Moberg, G. P. Biological response to stress: implications for animal welfare in The Biology of Animal Stress: Basic Principles and Implications for Animal Welfare (eds Moberg, G. P. & Mench, J. A.) 1–21 (CABI Publishing, 2000).18.Romero, M. L. & Butler, L. K. Endocrinology of stress. Int. J. Comp. Psychol. 20, 89–95 (2007).
Google Scholar 
19.Shutt, K., Setchell, J. M. & Heistermann, M. Non-invasive monitoring of physiological stress in the Western lowland gorilla (Gorilla gorilla gorilla): Validation of a fecal glucocorticoid assay and methods for practical application in the field. Gen. Comp. Endocrinol. 179, 167–177 (2012).CAS 
PubMed 
Article 

Google Scholar 
20.Heistermann, M. Non-invasive monitoring of endocrine status in laboratory primates: Methods, guidelines and applications. Adv. Sci. Res. 5, 1–9 (2010).Article 

Google Scholar 
21.Murray, C. M., Heintz, M. R., Lonsdorf, E. V., Parr, L. A. & Santymire, R. M. Validation of a field technique and characterization of fecal glucocorticoid metabolite analysis in wild chimpanzees (Pan troglodytes). Am. J. Primatol. 75, 57–64 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
22.Schwarzenberger, F. The many uses of non-invasive faecal steroid monitoring in zoo and wildlife species. Int. Zoo. Yearb. 41, 52–74. https://doi.org/10.1111/j.1748-1090.2007.00017.x (2007).Article 

Google Scholar 
23.Touma, C. & Palme, R. Measuring fecal glucocorticoid metabolites in mammals and birds: The importance of validation. Ann. N. Y. Acad. Sci. 1046, 54–74. https://doi.org/10.1196/annals.1343.006 (2005).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
24.Kersey, D. C. & Dehnhard, M. The use of noninvasive and minimally invasive methods in endocrinology for threatened mammalian species conservation. Gen. Comp. Endocrinol. 203, 296–306. https://doi.org/10.1016/j.ygcen.2014.04.022 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
25.Palme, R., Rettenbacher, S., Touma, C., El-Bahr, S. M. & Möstl, E. Stress hormones in mammals and birds: Comparative aspects regarding metabolism, excretion, and noninvasive measurement in fecal samples. Ann. N. Y. Acad. Sci. 1040, 162–171 (2005).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
26.Teskey-Gerstl, A., Bamberg, E., Steineck, T. & Palme, R. Excretion of corticosteroids in urine and faeces of hares (Lepus europaeus). J. Comp. Physiol. B. 170, 163–168 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
27.Millspaugh, J. J. & Washburn, B. E. Use of fecal glucocorticoid metabolite measures in conservation biology research: Considerations for application and interpretation. Gen. Comp. Endocrinol. 138, 189–199 (2004).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
28.Bahr, N. I., Palme, R., Möhle, U., Hodges, J. K. & Heistermann, M. Comparative aspects of the metabolism and excretion of cortisol in three individual nonhuman primates. Gen. Comp. Endocrinol. 117, 427–438. https://doi.org/10.1006/gcen.1999.7431 (2000).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
29.Whitten, P. L., Brockman, D. K. & Stavisky, R. C. Recent advances in noninvasive techniques to monitor hormone-behavior interactions. Am. J. Phys. Anthropol. 107, 1–23 (1998).Article 

Google Scholar 
30.Heistermann, M., Palme, R. & Ganswindt, A. Comparison of different enzymeimmunoassays for assessment of adrenocortical activity in primates based on fecal analysis. Am. J. Primatol. 68, 257–273 (2006).CAS 
PubMed 
Article 

Google Scholar 
31.Wheeler, B. C., Tiddi, B., Kalbitzer, U., Visalberghi, E. & Heistermann, M. Methodological considerations in the analysis of fecal glucocorticoid metabolites in tufted capuchins (Cebus apella). Int. J. Primatol. 34, 879–898 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
32.Pei, J.-C. et al. Disease surveillance, conservation, and management strategies for Taiwanese macaque (Macaca cyclopis) at Shoushan National Nature Park. 136 (Construction and Planning Agency. Ministry of the Interior, 2015).33.Hsu, M. J., Kao, C. C. & Agoramoorthy, G. Interactions between visitors and Formosan macaques (Macaca cyclopis) at Shou-Shan Nature Park, Taiwan. Am. J. Primatol. 71, 214–222 (2009).PubMed 
Article 

Google Scholar 
34.Lee, L. L., Wu, H. Y., Chang, S. W., Minna, J. H. & Chakraborty, C. Survey of Current Status of Taiwan Macaques 1–27 (Council of Agriculture, Executive Yuan, 2001).35.Pei, K. C. J., Chen, C. C., Lin, C. N. & Ju, Y. T. The study of population dynamic and health status of Taiwanese macaques in Shanshan National Nature Park., 175 (Shoushan National Nature Park, Construction and Planning Agency, Ministry of the Interior, 2016).36.Bellanca, R. U. et al. A simple alopecia scoring system for use in colony management of laboratory-housed primates. J. Med. Primatol. 43, 153–161 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
37.Whiting, D. A. Histology of the human hair follicle in Hair Growth and Disorders (eds Blume-Peytavi, U. et al.) 107–123 (Springer, 2008).38.Schneider, C. A., Rasband, W. S. & Eliceiri, K. W. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671–675 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Horenstein, M. G. & Bacheler, C. J. Follicular density and ratios in scarring and nonscarring alopecia. Am. J. Dermatopathol. 35, 818–826 (2013).PubMed 
Article 

Google Scholar 
40.Rangel-Negrín, A., Flores-Escobar, E., Chavira, R., Canales-Espinosa, D. & Dias, P. A. D. Physiological and analytical validations of fecal steroid hormone measures in black howler monkeys. Primates 55, 459–465 (2014).PubMed 
Article 

Google Scholar 
41.Pineda-Galindo, E., Cerda-Molina, A. L., Mayagoitia-Novales, L. & Matamoros-Trejo, G. Biological validations of fecal glucocorticoid, testosterone, and progesterone metabolite measurements in captive stumptail macaques (Macaca arctoides). Int. J. Primatol. 38, 985–1001 (2017).Article 

Google Scholar 
42.Palme, R. & Möstl, E. Measurement of cortisol metabolites in faeces of sheep as a parameter of cortisol concentration in blood. Int. J. Mammal. Biol. 62, 192–197 (1997).
Google Scholar 
43.Braga Goncalves, I. et al. Validation of a fecal glucocorticoid assay to assess adrenocortical activity in meerkats using physiological and biological stimuli. PLoS ONE 11, e0153161. https://doi.org/10.1371/journal.pone.0153161 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
44.Pei, K. J.-C., Lin, C. N. & Chen, C. C. Disease surveillance, conservation, and management strategies for Taiwanese macaque (Macaca cyclopis) at Shoushan National Nature Park. 133 (Construction and Planning Agency, Ministry of the Interior, R2015).45.Hsu, M. J. & Lin, J.-F. Troop size and structure in free-ranging Formosan macaques (Macaca cyclopis) at Mt. Longevity, Taiwan. Zool. Stud. Taipei 40, 49–60 (2001).
Google Scholar 
46.Graham, M. H. Confronting multicollinearity in ecological multiple regression. Ecology 84, 2809–2815 (2003).Article 

Google Scholar 
47.Allison, P. D. Multiple Regression: A Primer (Pine Forge Press, 1999).
Google Scholar 
48.Dohoo, I., Martin, W. & Stryhn, H. Model-building strategies in Veterinary Epidemiologic Research 553–578 (VER Inc, 2009).49.Bolker, B. et al. Generalized linear mixed models: A practical guide for ecology and evolution. Trends Ecol. Evol. 24, 127–135 (2009).PubMed 
Article 

Google Scholar 
50.Ganswindt, A., Palme, R., Heistermann, M., Borragan, S. & Hodges, J. Non-invasive assessment of adrenocortical function in the male African elephant (Loxodonta africana) and its relation to musth. Gen. Comp. Endocrinol. 134, 156–166 (2003).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Rimbach, R., Heymann, E. W., Link, A. & Heistermann, M. Validation of an enzyme immunoassay for assessing adrenocortical activity and evaluation of factors that affect levels of fecal glucocorticoid metabolites in two New World primates. Gen. Comp. Endocrinol. 191, 13–23 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.Novak, M. A. et al. Assessing significant ( > 30%) alopecia as a possible biomarker for stress in captive rhesus monkeys (Macaca mulatta). Am. J. Primatol. 79, 1–8. https://doi.org/10.1002/ajp.22547 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
53.Bernardez, C., Molina-Ruiz, A. & Requena, L. Histologic features of alopecias–part I: Nonscarring alopecias. Actas Dermosifiliogr. 106, 158–167. https://doi.org/10.1016/j.adengl.2015.01.001 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
54.Luchins, K. R. et al. Application of the diagnostic evaluation for alopecia in traditional veterinary species to laboratory rhesus macaques (Macaca mulatta). J. Am. Assoc. Lab. Anim. Sci. 50, 926–938 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
55.Werner, B. & Mulinari-Brenner, F. Clinical and histological challenge in the differential diagnosis of diffuse alopecia: Female androgenetic alopecia, telogen effluvium and alopecia areata-part II. An. Bras. Dermatol. 87, 884–890 (2012).PubMed 
PubMed Central 
Article 

Google Scholar 
56.Liyanage, D. & Sinclair, R. Telogen effluvium. Cosmetics 3, 13 (2016).Article 
CAS 

Google Scholar 
57.Alotaibi, M. K. Telogen effluvium: A review. Int. J. Med. Dev. Cties. 3, 797–801. https://doi.org/10.7759/cureus.8320 (2019).Article 

Google Scholar 
58.Arck, P. C. et al. Stress inhibits hair growth in mice by induction of premature catagen development and deleterious perifollicular inflammatory events via neuropeptide substance P-dependent pathways. Am. J. Pathol. 162, 803–814. https://doi.org/10.1016/S0002-9440(10)63877-1athology (2003).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
59.Horenstein, V.D.-P., Williams, L. E., Brady, A. R., Abee, C. R. & Horenstein, M. G. Age-related diffuse chronic telogen effluvium-type alopecia in female squirrel monkeys (Saimiri boliviensis boliviensis). Comp. Med. 55, 169–174 (2005).PubMed 
PubMed Central 

Google Scholar 
60.Coleman, K. et al. The correlation between alopecia and temperament in rhesus macaques (Macaca mulatta) at four primate facilities. Am. J. Primatol. 79, e22504. https://doi.org/10.1002/ajp.22504 (2017).Article 

Google Scholar 
61.Lutz, C. K. et al. Factors influencing alopecia and hair cortisol in rhesus macaques (Macaca mulatta). J. Med. Primatol. 45, 180–188. https://doi.org/10.1111/jmp.12220 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
62.Palme, R. Non-invasive measurement of glucocorticoids: Advances and problems. Physiol. Behav. 199, 229–243. https://doi.org/10.1016/j.physbeh.2018.11.021 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
63.Hoffman, C. L. et al. Immune function and HPA axis activity in free-ranging rhesus macaques. Physiol. Behav. 104, 507–514. https://doi.org/10.1016/j.physbeh.2011.05.021 (2011).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
64.Marty, P. R., Hodges, K., Heistermann, M., Agil, M. & Engelhardt, A. Is social dispersal stressful? A study in male crested macaques (Macaca nigra). Horm. Behav. 87, 62–68. https://doi.org/10.1016/j.yhbeh.2016.10.018 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
65.Takeshita, R. S. C., Bercovitch, F. B., Kinoshita, K. & Huffman, M. A. Beneficial effect of hot spring bathing on stress levels in Japanese macaques. Primates 59, 215–225. https://doi.org/10.1007/s10329-018-0655-x (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
66.Sheriff, M. J., Dantzer, B., Delehanty, B., Palme, R. & Boonstra, R. Measuring stress in wildlife: Techniques for quantifying glucocorticoids. Oecologia 166, 869–887 (2011).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
67.Cheng, H. C. et al. The Red List of Terrestrial Mammals of Taiwan, 2017. 35 (Endemic Species Research Institute, 2017).68.Chang, A.-M., Chen, C.-C. & Huffman, M. A. Entamoeba spp in wild formosan rock macaques (Macaca cyclopis) in an area with frequent human-macaque contact. J. Wildl. Dis. 55, 608–618 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar  More

In this multi-scale assessment of temporal variation in bird-window collisions, our predictions related to the diel timing of collisions were only partly supported. We predicted more casualties would occur during morning than other times of day, which should have resulted in our detection of the greatest number of fatal collisions on midday surveys and more non-fatal collisions during morning surveys than midday and evening surveys. However, greatest numbers of both fatal and non-fatal collisions were observed on morning surveys, indicating that more collisions occurred during overnight and early morning periods than mid-to-late morning and afternoon. As predicted, this diel pattern was consistent across seasons. Our predictions about monthly and seasonal patterns were also only partly supported. Unexpectedly, total collision mortality was highest in the spring migration month of May, and avian residency status interacted with season such that roughly equivalent high numbers of resident and migrant collisions occurred in spring, more resident than migrant collisions occurred in summer (and overall from Apr to Oct), and more migrant than resident collisions occurred in fall. Further, unlike many past studies, collision mortality of migrants was roughly equivalent in spring and fall migration.Diel collision patternsWe observed more fatal and non-fatal collisions in the morning than in midday and evening surveys combined, even though we included fewer morning surveys in diel analyses. These differences are likely conservative in that an even greater proportion of fatal collisions than we observed likely occurred overnight and in the early morning, but went undetected because of bias associated with observer detection and scavenger removal of carcasses. Concurrent work in our study area20 found that relatively inexperienced volunteers detected a slightly lower proportion of carcasses (0.69) than experienced surveyors (0.76). Because roughly 10% of morning surveys were conducted by less-experienced volunteers and all midday or evening surveys were conducted by full-time technicians or the authors, we likely missed more carcasses on morning surveys. Additionally, most scavenging events (68%) were at night20, so bird carcasses from overnight collisions were the least likely to persist until the subsequent survey. Thus, underestimation of fatal collisions in the preceding interval was almost certainly greater for morning surveys than for midday and evening surveys.A prevailing hypothesis for why daytime bird-window collisions occur is that birds making local (e.g., foraging) movements fail to perceive a barrier when flying toward objects either on the other side of glass or reflected on a glassy surface8,15. Under this hypothesis, daytime collisions for both residents and migrating birds at stopover locations would be expected to occur most frequently when birds are most active, which is typically near dawn regardless of season. Our finding of the greatest number of collision fatalities on morning surveys circumstantially supports the above hypothesis and expectation, as do past descriptive studies of diel variation in bird-window collisions26,27,38,39. However, our study design did not allow differentiation between nocturnal and early morning collisions, and a nontrivial proportion of carcasses detected on morning surveys likely represented collisions from the preceding overnight period. Nighttime collisions may occur at any structural component not easily detectable at night (i.e., they are not limited to glass surfaces), and can be exacerbated by artificial light emission that attracts and disorients migrating birds36,37,44,45. Nonetheless, the observation of more non-fatal collisions (including directly witnessed collisions) during morning than midday or evening surveys does strongly suggest that morning carcass counts included many collisions that occurred near or shortly after dawn.A potential limitation of our study regarding time-of-day analyses is the longer interval between evening and morning surveys than between other survey periods. If collisions occurred uniformly or randomly in time, we would expect to find more bird carcasses during morning surveys due to the longer preceding time interval. However, as described above, the early morning peak observed for non-fatal collisions (Fig. 1), which are less persistent than carcasses and therefore do not accumulate over time, suggests that collisions do not have a uniform or random temporal distribution and that a real peak in collision frequency occurs sometime shortly prior to when we conducted morning surveys. Another possible bias is that we conducted surveys during fixed time periods rather than adjusting them to seasonally varying sunrise and sunset times. This limitation would have most strongly affected our summer results, potentially inflating summer morning counts as a result of collisions for both early morning and late evening (the two periods when birds are most active) being grouped into the same survey period. However, this limitation was unlikely to greatly influence our conclusions about diel collision patterns because: (1) these patterns were fairly consistent across seasons, suggesting a relatively small influence of varying sunrise and sunset times, and (2) sunrise and sunset times vary by only a few hours over the course of the entire year whereas the broad time periods for which we compared collisions (overnight, morning, afternoon) consist of longer lengths of time. Further research is needed to identify the exact timing of collisions, including during overnight periods, and this could be accomplished with collision surveys conducted more frequently during the day and night or remote detection methods, such as video cameras, motion-triggered still cameras, microphones that record sounds of impact, and glass-mounted pressure sensors that detect vibrations from collision impact.We predicted diel collision patterns would not vary seasonally because the pattern of birds being most active near sunrise and sunset is relatively consistent across seasons. Hager and Cosentino46 provide excellent guidelines for conducting bird-window collision surveys, but their recommendation to conduct surveys in mid-to-late afternoon is based on summer monitoring that found mortality to peak between late morning and early afternoon in Illinois, USA29. We suspect differences in diel patterns between that study and ours relate to geographic variation and/or seasonal sampling coverage, as our larger sample of surveys included spring and fall migration in addition to summer. Although many collision-prone species migrate nocturnally, the diel collision peak for migrants could still occur in early morning because nocturnally-migrating birds often set-down into stopover habitats during early morning47,48 and may be most susceptible to collisions at this time. There could be subtle seasonal variation in diel collision patterns that we failed to detect; however, the majority of collisions across seasons appear to occur near or before dawn (see also26,27,38,39). In combination with the previous study showing that scavenging peaks overnight20, we suggest that conducting daily collision surveys in the morning could result in the least biased mortality estimates, especially in urbanized areas where humans (e.g., cleaning crews) may remove carcasses in the early morning. As noted by Hager and Cosentino46, further research is needed to identify how the optimal survey time is influenced by factors such as geography, the bird community, animal scavengers, and removal of bird carcasses by humans.Monthly and seasonal collision patternsWe expected more collisions in fall than other seasons because bird populations in North America are larger after summer breeding and include higher proportions of young birds that have less experience with flight, migration, and human structures. Also, numerous studies have found the greatest window collision mortality in fall, a pattern driven largely by migrant birds 12,18,23,27,31,40,49. Contrary to expectation, both raw counts and bias-adjusted estimates of collision fatalities were highest in the spring migration month of May and higher overall in spring than fall. This pattern resulted from a high number of both migrant and resident bird collisions. In fact, we found that roughly equal numbers of resident and migrant collisions occurred in spring. When considering migrating individuals only, we found roughly similar numbers of collisions in spring and fall migration, a finding that also is unexpected in light of past studies. Notably, two other studies that found that a large proportion of total collisions consisted of resident birds30,50 also documented a seasonal pattern of total collisions that was less skewed toward fall. An explanation for the relatively large amount of total spring mortality, and for our finding that migrant mortality was roughly similar in spring and fall, is that some long-distance migrants follow elliptical migration paths where migration routes in fall are farther east than in spring51,52, such that in central North America, numbers of some species are greatest during spring migration. This explanation is supported by our observation of some elliptical migrants colliding during spring but not fall (e.g., Swainson’s Thrush [Catharus ustulatus]). Our study adds further nuance to the understanding of seasonal variation in bird collisions and exemplifies the need to study bird-window collisions in a wider array of geographic contexts to allow region-relevant management recommendations.Our predictions regarding avian residency status were only partly supported; more migrants than non-migrants were indeed killed during fall migration. In spring, however, roughly equal numbers of each group collided, and far more residents than migrants collided in summer (and overall from Apr–Oct). This result does not account for varying abundance (overall and seasonally) of each species group, so it does not necessarily imply that resident species are more vulnerable to, or at greater risk of, window collisions relative to their abundance or period of residency in our study area. However, the finding of a high proportion of non-migrant casualties was still unexpected given that previous similar studies have almost universally reported higher mortality among migrants15,25,26,39,49,53,54, although most sampled during migratory periods only. Even with our individual-based residency designations, we may have slightly underestimated migrant mortality because all individuals of some migratory species were classified as unknown due to overlapping resident and migratory periods. However, even if all unknown individuals were migrating, there were too few birds in this category (22 of 341 [7%] total carcasses) to change our conclusions regarding the migrant-resident comparisons. Anecdotally, many spring and summer collision fatalities were recently fledged juveniles (n = 26 [25%] in May; n = 17 [30%] in June), clearly indicating that some collision victims were indeed not migrating, and therefore, that the high number of resident collisions is not an artifact of our classification system. Moreover, we did not observe collisions of any migrant individuals during June, even though a few species migrate through our study area in small numbers during this period (e.g., shorebirds [order Charadriiformes] and some tyrant flycatchers [family Tyrannidae])55. It is possible, however, that some resident individuals were undergoing post-breeding dispersal at the time of collision, as evidenced by a small late-June peak of Tufted Titmouse (Baeolophus bicolor) and Carolina Chickadee (Poecile carolinensis) collision victims with brood patches (TJO unpublished data).Although our sampling captured the peak months of spring and fall migration in our study region, we undoubtedly missed some early-spring migrants before April and late-fall migrants after October. However, greater than 10 years of near-daily collision surveys at one of the most collision-prone buildings in our study (TJO unpublished data) suggests this number of missed collisions was relatively small. Specifically, total collisions at this building (including residents and migrants) dropped from an average of  > 8 birds in October to  More

Global Health Program, Smithsonian Conservation Biology Institute, Washington DC, USAJames M. Hassell & Dawn ZimmermanDepartment of Epidemiology of Microbial Disease, Yale School of Public Health, New Haven, CT, USAJames M. Hassell & Dawn ZimmermanCentre for Biodiversity & Environment Research (CBER), Department of Genetics, Evolution and Environment, University College London, London, UKTim Newbold & Lydia H. V. FranklinosDepartment of Ecology & Evolutionary Biology, Princeton University, Princeton, NJ, USAAndrew P. DobsonSanta Fe Institute, Santa Fe, NM, USAAndrew P. DobsonWalter Reed Biosystematics Unit (WRBU), Smithsonian Institution Museum Support Center, Suitland, MD, USAYvonne-Marie LintonDepartment of Entomology, Smithsonian National Museum of Natural History, Washington DC, USAYvonne-Marie LintonWalter Reed Army Institute of Research (WRAIR), Silver Spring, MD, USAYvonne-Marie LintonMarine Disease Ecology Laboratory, Smithsonian Environmental Research Center, Edgewater, MD, USAKatrina M. Pagenkopp Lohan More

Plant growthPrimary root length and rosette size were estimated. Control root length, measured until 21 DAG, was lightly minor in ddc vs WT (Fig. 1A). Cd differentially inhibited root growth in the two samples: at 21 DAG, 25 and 50 µM Cd-treated roots were 1.2 and 2.2 fold shorter than control roots in ddc, while in the WT Cd-treated roots were 1.8 and 2.8 fold shorter than control ones (Fig. 1A). Consequently, at 21 DAG root of Cd-treated samples was longer in ddc vs WT, particularly at the lowest Cd concentration.Figure 1(A) Primary root length (B) Picture of rosette leaf series and (C) rosette leaf area (cm2) of WT and ddc plants of A. thaliana, germinated and grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl). Root length was monitored up to 21 days after germination (DAG) every two days from germination. The results represent the mean value (± SD) of three independent biological replicates (n = 45). Asterisks indicate significant pairwise differences using Student’s t-test (*P ≤ 0.05; ** P ≤ 0.01; *** P ≤ 0.001), performed between ddc vs WT subjected to the same treatment. Bars, 0.5 cm.Full size imageRosette size was estimated at 21 DAG, corresponding to the period necessary for its full development11, by evaluating leaf number and area. Control plants of both ddc and WT exhibited a complete leaf series, although most leaves resulted smaller in ddc (Fig. 1B,C). Cd affected rosette development reducing leaf number and area, less in ddc than WT, resulting into a higher leaf area and/or number in ddc under both Cd concentrations (Fig. 1B,C).Gene expression profileRNA-Seq analysis provided an overview of gene expression profile of Cd-treated and control plants of both ddc and WT. The following comparisons were performed: ddc vs WT under control (Ctrl) conditions (ddc vs WT-Ctrl) and 25 and 50 µM Cd treatment (ddc vs WT-25 µM Cd; ddc vs WT-50 µM Cd); 25/50 µM Cd-treated vs Ctrl in ddc (25 µM Cd vs Ctrl-ddc; 50 µM Cd vs Ctrl-ddc); 25/50 µM Cd-treated vs Ctrl in the WT (25 µM Cd vs Ctrl-WT; 50 µM Cd vs Ctrl-WT).After DEGs identification (see Supplementary Fig. S1 online) 14 of them were analysed through qRT-PCR to validate transcriptomic analysis (see Supplementary Fig. S2 online). Results were fully consistent with RNA-seq data. Gene Enrichment analysis was also performed, evidencing that Cd strongly impacted on transcriptome in both ddc and WT, but in a largely different way (see Supplementary Figs. S3–S9 online). Notwithstanding, a common aspect was that in both ddc and WT the genetic pathways (GPs) more impacted by Cd dealt with photosynthesis, stress responses and hormone biosynthesis and signalling.Expression pattern of genetic pathways related to hormonesIn view of hormones pivotal role in plant development and stress response and considering the assessed epigenetic control on their action and signalling12, in this work we analysed in depth how the expression pattern (EP) of hormone-related GPs was modulated in ddc vs WT under Cd stress. The most relevant differences are discussed.AuxinsUnder control conditions, GPs related to auxin biosynthesis showed comparable EP in ddc and WT and no DEGs were detected (Fig. 2A). 25 µM Cd induced significant changes only in ddc resulting into: i) TAA1 and YUC5 downregulation along indole-3-pyruvic acid (IPA) pathway; ii) CYP71A13 and NIT2 overexpression along indole-3-acetaldoxime (IAOX) auxiliary pathway, while CYP79B3 was downexpressed (Fig. 2A). Differently, 50 µM Cd induced similar changes in ddc and WT consisting in: (i) a downexpression of YUC2 along IPA pathway in both samples and YUC5 and YUC9 in ddc and WT, respectively; (ii) overexpression of CYP71A12, CYP71A13, NIT2, NIT4 and downexpression of CYP71A16 along IAOX pathway in ddc and WT (Fig. 2A).Figure 2Genes differentially expressed (DEGs) along the pathway of (A) auxin biosynthesis, auxin conjugation, (B) indole-3-acetyl-amino acid biosynthesis, (C) methyl-indole-3-acetate interconversion and (D) auxin signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageAuxin level and homeostasis also depend on its oxidative degradation, conjugation and methylation13. Under control conditions, GPs related to auxin conjugation and methylation showed comparable EPs in ddc and WT and no DEGs were detected (Fig. 2B,C), but were differentially impacted by Cd, mainly at 25 µM concentration. Namely, at 25 µM Cd several genes related to auxin conjugation (GH3.3, GH3.17, YDK1) and methylation (MES7, MES17) were downregulated in ddc, whereas in WT only MES18, involved in methyl-indole-3-acetate production, was downregulated (Fig. 2B,C). At 25 µM Cd most of the above genes were downregulated in ddc vs WT (Fig. 2B,C), while at 50 µM Cd only two genes, working in auxin methylation were differentially modulated in ddc and WT (Fig. 2B,C).Under control conditions, GP related to auxin signalling exhibited similar EP in ddc and WT (Fig. 2D). Differences were induced by Cd. Namely, at 25 µM Cd, AUX/IAA family genes, which acts in signalling repression, were globally downregulated in ddc, while in WT the repressor IAA34 was overexpressed (Fig. 2D). Differently, 50 µM Cd effects on ddc and WT were quite similar, dealing with AUX/IAA family genes downregulation (Fig. 2D). In ddc vs WT comparisons, at 25 µM Cd above genes were downregulated, while no differences were found at 50 µM Cd (Fig. 2D). Moreover, and somehow unexpectedly, following 50 µM Cd it was observed in both ddc and WT a downregulation of several SAURs members, belonging to a large family of auxin responsive genes, which in turn can also have an impact on auxin pathway14 (Fig. 2D). Interestingly, such effect was more pronounced in ddc than WT. However, it must be mentioned that, although most of them are induced by auxin, several other hormones and co-factors acts upstream SAUR genes, regulating their activity in response to both endogenous stimuli and environmental cues14.In summary, in both ddc and WT, Cd induced: i) a downregulation of IPA pathway, which is the main auxin biosynthetic pathway15 and a simultaneous upregulation of IAOX auxiliary biosynthetic pathway; ii) an enhancement of hormone signalling. However, in the WT such effects occurred only at 50 μM Cd. Moreover, in ddc Cd also induced a downregulation of GPs related to auxin conjugation.CytokininsIn all comparisons, GPs related to CKs biosynthesis showed similar EPs, unless for the downregulation in 50 µM Cd-treated WT vs Ctrl of IPT5, encoding rate-limiting enzyme along the pathway16 (Fig. 3A).Figure 3Genes differentially expressed (DEGs) along the pathway of (A) trans-zeatin biosynthesis, (B) CKs degradation, (C) CKs N7- and N9-glucoside biosynthesis, (D) CKs O-glycosylation and (E) CKs signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageMajor differences were observed for GPs related to CKs catabolism and conjugation, occurring through cleavage by oxidation and glycosylation, respectively16. Under control conditions, these GPs also exhibited similar EPs in ddc and WT (Fig. 3B-D). At both 25 and 50 µM Cd a downregulation of CKX5 and CKX6, encoding cytokinin-oxidases, occurred only in ddc (Fig. 3B). Differently, Cd impact on GPs related to CKs N-glycosylation was almost comparable in ddc and WT, resulting into the overexpression of two different genes working in N7- and N9-glycosylation pathways at 25 µM Cd, and one gene at the higher concentration (Fig. 3C).Note that at 25 µM Cd the above genes were both overexpressed in ddc vs WT, while no differences were found at 50 µM Cd (Fig. 3C). Cd impact on GP related to cytokinin O-glycosylation was major, especially in ddc, involving at 25 µM Cd the overexpression of seven genes along this pathway compared to Ctrl (Fig. 3D), and one gene in WT 25 µM Cd vs Ctrl (Fig. 3D). By contrast, at 50 μM Cd the expression pattern along this pathway was similar in ddc and WT, being characterised by the upregulation of the same seven genes above mentioned and the downregulation of AT5G38010 (Fig. 3D). Finally, at 25 μM Cd, five genes along these pathways resulted upregulated in ddc vs WT while a similar EP occurred in ddc and WT under 50 μM Cd treatment (Fig. 3D).Concerning GP involved in CKs signalling, under control conditions A-ARRs, encoding negative regulators of CKs signalling17, were downregulated and signalling was likely enhanced in ddc vs WT (Fig. 3E). 25 µM Cd induced a downregulation of ARR11 A-type ARRs and B-ARR family ARR10 transcription factors, which control primary plant response to CKs, only in ddc (Fig. 3E). Whereas, at 50 µM Cd both ddc and WT showed A-ARR downregulation, supposedly leading to pathway upregulation (Fig. 3E). At 25 µM Cd AHP1, encoding positive regulators of CKs signalling18, was downregulated in ddc vs WT, while ARR17 was overexpressed, suggesting that signalling was downregulated also in ddc vs WT (Fig. 3E). No differences occurred between ddc and WT at 50 µM Cd (Fig. 3E).In summary, transcriptomic analysis evidenced that GP related to the biosynthesis of trans-zeatin, the most relevant CK, was negatively affected by Cd only in the WT at 50 μM Cd. In response to Cd, GPs related to CKs inactivation were enhanced in both ddc and WT, but in ddc a downregulation of GP related to CKs cleavage also occurred. Finally, hormone signalling was differentially modulated by Cd in relation to both the sample (ddc vs WT) and heavy metal concentration, resulting into a downregulation at 25 µM Cd only in ddc and an enhancement in ddc and WT at 50 µM Cd.GibberellinsUnder control conditions, GPs related GAs biosynthesis showed similar EPs in ddc vs WT (Fig. 4A). 25 μM Cd induced in ddc: i) a downregulation of GA2 encoding the ent-kaurene synthase, a pivotal enzyme along the early GAs biosynthetic pathways to synthetize GA12; ii) a downregulation of GA4, a key gene of GAs biosynthesis, along which bioactive GAs are synthetized19 (Fig. 4A). No Cd-induced modulation was observed in the WT (Fig. 4A). On the contrary, at 50 μM Cd both ddc and WT showed a downregulation of GA5 (Fig. 4A). Finally, in ddc vs WT the only difference dealt with GA2 downregulation at 25 µM Cd (Fig. 4A).Figure 4Genes differentially expressed (DEGs) along the pathway of (A) GAs biosynthesis, (B) GAs inactivation and (C) GAs signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageGPs controlling GAs inactivation also showed a comparable transcriptional pattern in ddc and WT under control conditions, and no DEGs were detected (Fig. 4B). 25 μM Cd induced a downregulation of DAO2 and AOP1, encoding GA2ox enzymes, only in ddc (Fig. 4B). Accordingly, in ddc vs WT these genes were downregulated only at the lowest Cd concentration (Fig. 4B).Under control conditions, also GP related to GAs signalling was not differentially modulated in ddc vs WT (Fig. 4C). A downregulation of genes encoding DELLA proteins, which act as repressors20, was induced only in ddc by 25 µM Cd (Fig. 4C) and in both ddc and WT at 50 µM Cd (Fig. 4C), suggesting an enhancement of hormone signalling. Consequently, DELLA-codifying genes resulted downregulated also in ddc vs WT only at the lowest Cd concentration (Fig. 4C).In summary, in ddc the lowest Cd treatment negatively affected GPs related to GAs biosynthesis but, at the same time, hormone signalling resulted enhanced. In the WT similar effects were observed only at 50 µM Cd.Jasmonic acidUnder control conditions, six genes along the GP related to JA biosynthesis were downregulated in ddc vs WT (Fig. 5A). 25 µM Cd induced in ddc a downregulation of this GP, except for LOX4 upregulation (Fig. 5A) and a downregulation involving eight genes in WT (Fig. 5A). At 50 μM, Cd effects were limited to LOX5 downregulation and LOX4 and OPR1 upregulation in ddc and LOX5 and AOS downregulation in WT (Fig. 5A). No Cd–induced differences were found in ddc vs WT (Fig. 5A).Figure 5Genes differentially expressed (DEGs) along the pathway of (A) JA biosynthesis, (B) JA signalling, (C) ABA biosynthesis, (D) ABA degradation, (E) ABA glucose ester biosynthesis and (F) ABA signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageConcerning the JA signalling-related GP, in control conditions JAZ5 gene, encoding a protein acting as repressor21, was downregulated in ddc vs WT, highlighting a signalling enhancement (Fig. 5B). Interestingly, JAZ10 and JAZ9 were differentially impacted by 25 µM Cd in ddc and WT resulting upregulated and downregulated, respectively. (Fig. 5B). At 50 µM Cd, JAZs were overexpressed in both ddc and WT (Fig. 5B). No differences were detected in ddc vs WT exposed to Cd (Fig. 5B).Globally, Cd negatively impacted on the GP related to JA biosynthesis especially in WT. Under Cd treatment hormone signalling was downregulated more in ddc than in WT, whatever concentration was applied.Abscisic acidUnder control conditions, ABA biosynthesis-related GP showed comparable EP in ddc and WT, and no DEGs were detected (Fig. 5C). The only significant Cd effect dealt with NCED3 downregulation both in ddc and WT, regardless of applied concentration (Fig. 5C). Under control conditions, also the GPs related to ABA catabolism showed a comparable EP in ddc and WT and no DEGs were detected (Fig. 5D), but Cd differentially impacted on CYP genes, involved in phaseic acid degradative production22. Namely, at 25 μM Cd, CYP707A3 was downregulated only in ddc (Fig. 5D). Moreover, also CYP707A2 appeared downregulated in ddc vs WT (Fig. 5D). At 50 µM Cd it was observed an upregulation of both CYP707A2 and CYP707A4 in ddc, and of only CYP707A4 in WT (Fig. 5D).Under control conditions, GP related to ABA inactivation through glucose conjugation showed similar EP in ddc vs WT (Fig. 5E). 25 μM Cd determined AT4G15260 upregulation and UGT71C3 downregulation only in ddc (Fig. 5E). At 50 μM, Cd equally impacted on ddc and WT, resulting into UGT71C1 and UGT2 downregulation, AT5G49690, UGT71B5, UGT71B6 upregulation and, limited to WT, AT4G15260 upregulation (Fig. 5E). No differences were highlighted in ddc vs WT exposed to Cd (Fig. 5E).Under control conditions, the GP related to hormone signalling also presented a comparable EP in ddc and WT (Fig. 5F). In ddc, 25 µM Cd impact on this GP appeared rather complex, resulting in an upregulation of PYL3, encoding ABA receptor, and a downregulation of PP2Cs (PP2CA and HAI1) encoding negative regulators of ABA signalling23. Moreover, ABI5, codifying a key transcription factor in ABA signalling24 belonging to AREBs/ABFs family, was upregulated. However, SnRK2.7 gene, codifying a protein which activate the AREBs/ABFs transcription factors24, was downregulated. Based on the prominent role of SnRK2s in plant response to ABA, it is likely that at 25 µM Cd ABA signalling was downregulated in ddc (Fig. 5F). Instead, 25 µM Cd determined in WT the upregulation of PYL6 and the downregulation of PP2Cs, suggesting an enhancement of ABA signalling (Fig. 5F). At 50 µM Cd, PYL6 was upregulated and SnRK2.7 downregulated in both ddc and WT, while HAI1 was downregulated only in ddc (Fig. 5F). When comparing ddc vs WT, at 25 µM Cd only SnRK2.7 was downregulated, while no differences occurred at 50 μM Cd (Fig. 5F).In summary, Cd determined a slight downregulation of GP related to ABA biosynthesis in both ddc and WT regardless of its concentration. ABA catabolic pathway was lightly downregulated in ddc at 25 µM Cd but upregulated in both samples at 50 μM Cd. At the transcriptomic level, ABA signalling featured as enhanced in WT and downregulated in ddc regardless of Cd concentration.EthyleneAlong GP related to ethylene biosynthesis, under control conditions ACS8 and ACS11 were upregulated in ddc vs WT (Fig. 6A). 25 µM Cd determined ACS8 and ACO5 downregulation in ddc and ACS4 upregulation in the WT (Fig. 6A). 50 μM Cd induced ACS7 upregulation and ACO5 downregulation in both ddc and WT and ACS2 and ACS11 overexpression only in WT (Fig. 6A). Finally, the only Cd-induced difference in ddc vs WT dealt with ACO1 downregulation at 25 µM Cd (Fig. 6A).Figure 6Genes differentially expressed (DEGs) along the pathway of (A) ethylene biosynthesis, (B) ethylene signalling and (C) SA signalling in ddc and WT plants identified through a transcriptomic approach. For each comparison, the log2(fold change) of the analysed DEGs was shown in orange and in blue for the upregulated and downregulated genes, respectively. Plants were grown for 21 DAG in long day condition: (i) on growth medium added with 25 or 50 µM Cd; (ii) on growth medium without Cd as control (Ctrl).Full size imageGP related to ethylene signalling showed similar EP in ddc and WT both under control conditions (Fig. 6B) and at 25 µM Cd, except for the upregulation of ETR2, encoding ethylene receptor25 in WT (Fig. 6B). At 50 µM Cd, both ddc and WT exhibited ETR2 and ERF1 overexpression suggesting an upregulation of ethylene signalling (Fig. 6B); no differences occurred in ddc vs WT (Fig. 6B).In summary, in control conditions GP related to ethylene biosynthesis was upregulated in ddc vs WT. Cd determined a downregulation and upregulation of this GP in ddc and WT, respectively. Concerning hormone signalling, at the highest Cd concentration in both ddc and WT an upregulation of this GP occurred.Salicylic acidRegarding SA, only the GP related to signalling resulted differentially expressed. Under control conditions, the GP related to SA signalling showed similar EP in ddc and WT. However, PR1, a useful molecular marker for the systemic acquired resistance (SAR) in response to pathogens26, was downregulated in ddc vs WT (Fig. 6C). 25 μM Cd induced a downregulation of genes codifying TGA10 transcription factor only in the ddc (Fig. 6C). Whereas, 50 μM Cd induced a downregulation of PRB1 in both ddc and WT and of TGA8 only in ddc (Fig. 6C). In ddc vs WT, differences were found only at 25 µM Cd, with the downregulation of TGA10 and CAPE3 (Fig. 6C).Altogether, these results evidenced a Cd-induced downregulation of this GP, likely resulting in an impairment of hormone signalling in both WT and ddc, but in the latter this effect already occurred at the lowest Cd concentration.Phytohormone levelBased on the major effects induced by 25 µM Cd treatment, hormone quantification was carried out on plants exposed to this concentration, compared to untreated control plants.Under control conditions, IAA amount was higher in ddc than WT, although not significantly. After Cd treatment, a decreasing trend was observed only in WT, resulting into a significant lower level as compared to ddc (Fig. 7A).Figure 7(A) IAA, (B–G) CKs, (H–R) GAs, (S) JA, (T) ABA, (U) SA and (V) SAG amount in A. thaliana ddc mutant and WT plants grown in Ctrl conditions and treated with 25 µM Cd estimated by GC–MS. The results represent the mean value (± SD) of three independent biological replicates. Statistical analysis was performed by using two-way ANOVA with Tukey post hoc test (P ≤ 0.05) after Shapiro–Wilk normality test. Means with the same letter are not significantly different at P ≤ 0.05.Full size imageConcerning CKs, both biological active (tZ) and inactive conjugate (tZR, cZR, tZOG, cZOG, iPR) forms were analysed (Fig. 7B–G). Under control conditions, all analysed CKs were present in ddc, but CKs conjugate forms and above all O-glycosylated exhibited the highest levels (Fig. 7B–G). By contrast, in WT tZ was not detectable and all the other CKs forms exhibited a lower level compared to the mutant, which appeared significant for tZR and cZOG (Fig. 7B–G). Following Cd treatment, CKs levels increased in ddc, except for tZ decrease. In WT, the unique Cd effect dealt with tZR increase and tZOG decrease. Consequently, under Cd treatment the level of all CKs forms remained higher (from 0.25 to 3 times) in ddc than in WT (Fig. 7B–G).Concerning GAs, precursors (GA9, GA19, GA20), biologically active forms (GA1, GA3, GA4, GA7) and catabolites (GA8, GA34, GA29, GA51) were analysed (Fig. 7H–R). Under control conditions, GA19 amount was significantly higher in ddc vs WT, while the amount of GA20, the other in serie precursor of hydroxylated forms, was comparable between the samples (Fig. 7H,I). Following Cd treatment, GA19 amount significantly increased in WT while a slight downtrend occurred in ddc, leading to comparable values in the two samples. The same trend was observed, but at less extent, for GA20 (Fig. 7H,I). Consistently, under control conditions also the amount of the active hydroxylated forms GA1 and GA3 was higher in WT than in ddc (Fig. 7J,K). Following Cd treatment, a decrease of their amount was detected only in WT, globally leading to a higher level of these GAs in ddc mutant compared to the WT (Fig. 7J,K). In addition, in ddc mutant also the related catabolites GA8 and GA29 were globally lower than in WT, under both control conditions and Cd treatment (Fig. 7L,M).Differences were observed also for GA9, precursor of non-hydroxylated GAs: under Cd treatment its amount decreased in the WT and was instead induced in ddc mutant, resulting in a quite comparable value between the two samples (Fig. 7N). Consistently, the amount of active non-hydroxylated forms, GA4 and GA7, increased under Cd treatment only in ddc mutant; also in this case, at the end of heavy metal treatment, comparable values were detected in ddc and WT (Fig. 7O,P) In agreement with these results, following Cd treatment, the amount of catabolites GA51 and GA34 did not change in the WT, whereas in ddc it increased and decreased, respectively (Fig. 7Q,R).As evident in Fig. 7S–V, differences were reported also for JA, ABA, SA and its predominant inactive conjugate, SA 2-O-β-D-glucoside (SAG). Under control conditions both JA and ABA amount was significantly lower in ddc vs WT and significantly decreased following Cd treatment only in the WT (Fig. 7S,T). Notwithstanding, under such condition the ABA amount remained lower in ddc than in WT while JA values were comparable in the two samples due its light, but not significant, increase in ddc (Fig. 7S,T). By contrast, under control conditions both SA and SAG amounts were significantly higher in ddc than in WT (Fig. 7U,V). Following Cd treatment, their amounts significantly decreased more in ddc than in WT, leading to an opposite condition (Fig. 7U,V).Testing of the involvement of SUPPRESSOR OF DRM1 DRM2 CMT3 (SDC) gene in ddc response to CdFinally, we planned to inquire on the possible involvement of SDC gene in the response of ddc triple mutant to Cd exposure. Indeed, it has been reported that in ddc mutant the misexpression of such gene, which encodes a F-Box protein, is ultimately responsible of the developmental phenotypes of ddc, such as curled leaves and reduced growth, as evidenced by its reversion in the drm1 drm2 cmt3 sdc quadruple mutant27. Note that in the WT SDC is silenced, being methylated in all its sequence contexts because of the redundant action of DRM2 and CMT3 enzymes. By contrast, in ddc, where DRM2 and CMT3 expression is silenced, the loss of non-CG methylation in the promoter region of SDC F-box gene determines its overexpression27.According to the above mentioned data27, we firstly verified that under control conditions SDC resulted silent in the WT and overexpressed in ddc also in our transcriptomic analysis (data not shown, complete raw transcriptomic data are available at NCBI SRA under the BioProject accessionPRJNA641242;https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA641242). Moreover, our data also showed that at the transcriptomic level SDC is not modulated by Cd since its expression level did not significantly change in ddc nor in WT whatever heavy metal concentration was applied.Thereafter, we tested the involvement of SDC in the growth response of ddc mutant under Cd exposure, by monitoring primary root length of Arabidopsis thaliana WT, ddc and sdc plants grown under the following conditions: (i) on growth medium without Cd as control (Ctrl) (ii) on a medium supplemented with 25/50 µM Cd; (iii) limited to the WT and sdc mutant plants, on a medium supplemented with 25/50 µM Cd plus 15 µM 5-Azacytidine (5-Aza), an inhibitor of DNA methylation applied in order to mimic the hypomethylated state of ddc mutant.Under control conditions, all three samples showed a similar root length. However, at 21 DAG, root was lightly shorter in ddc vs WT, while sdc displayed an intermediate length (Fig. 8A). Under both Cd treatments, roots were averagely longer in ddc than in WT. Again, sdc roots exhibited an intermediate length, more similar to WT than ddc (Fig. 8 B,C). Interestingly, WT and sdc plants treated with Cd plus 5-Aza had longer roots than the plants treated only with Cd, and quite comparable to ddc roots exposed to Cd (Fig. 8 D,E).Figure 8Primary root length of WT, sdc and ddc plants of A. thaliana, germinated and grown in long day condition (A) on growth medium without Cd as control (Ctrl), (B,C) on a medium supplemented with 25/50 µM Cd, (D,E) limited to WT and sdc plants, on a medium supplemented with 25/50 µM Cd plus 15 µM 5-Azacytidine (5-Aza). Root length was monitored up to 21 days after germination (DAG) every two days from germination. The results represent the mean value (± SD) of three independent biological replicates (n = 45). Statistical analysis was performed between samples at the same growth stage, by using two-way ANOVA with Tukey post hoc test (P ≤ 0.05) after Shapiro–Wilk normality test. Means with the same letter are not significantly different at P ≤ 0.05.Full size image More

Genotyping control and datasets creationsAfter filtering the initial raw dataset of 1,071 goats for poor quality genotype and related animals we obtained the dataset used for the haplotype sharing analysis, which included 980 animals and 48,396 SNPs. This dataset was further pruned for linkage disequilibrium (LD; r2  More

MDE hypothesisGradient features of species diversity of plant communities refer to regular changes in species diversity along a gradient of environmental factors at the community level12,43. The altitudinal gradient includes gradient effects of multiple environmental factors. It is therefore important to study altitudinal patterns of species diversity to reveal changes in biodiversity along environmental gradients. The width and range of species distribution along geographical gradients reflect the ecological adaptability, diffusivity, and evolutionary history of species44. To some extent, geographical distribution patterns of species diversity can be interpreted as outcomes of synthetic actions across altitudinal gradients resulting from eurychoric species with a greater distribution width and stenochoric species with a smaller distribution range along geographical gradients44. Hence, the MDE, environmental gradient, distribution area, human disturbance, and habitat heterogeneity all have effects on the vertical distribution patterns of species diversity45,46. According to the MDE hypothesis, there is overlap in the distribution range of species along altitudinal gradients; the highest overlap intensity occurs at intermediate elevations. There is relatively low overlap intensity at low and high elevation areas, and the peak values of species diversity occur at intermediate elevations.In this study, forest ecosystems on the Loess Plateau were separated into three communities: tree, shrub, and herb communities. The altitudinal patterns and factors that influence species diversity of mountainous vegetation were then determined at the plant community level related to the form and family of the plant species. We discovered that the family numbers of the herb, shrub, and tree communities reached their greatest values at intermediate, intermediate, and lower elevations, respectively. We also discovered that correlations of species diversity indices with elevations conformed to unimodal change patterns for herb, shrub, and tree communities, which presented their greatest values at higher, lower, and intermediate elevations, respectively. This showed that MDE is another important factor that affects the distribution patterns of species diversity along regional altitudinal gradients, in addition to temperature, precipitation, and the terrain.A large number of studies have already verified that MDE is a significant mechanism that influences the gradient patterns of species diversity. MDE not only functions along altitudinal gradients but also acts along latitudinal and temporal gradients28,29,30,47. However, the effects of MDE on species diversity patterns are highly controversial. Some studies considered MDE to be the main factor that results in unimodal patterns of species diversity47,48, whereas other studies affirmed that the effects of MDE are smaller in contrast to the functions of the distribution area, environmental gradient, and other factors29,30. Besides MDE, other factors may also lead to unimodal vertical gradient patterns of species diversity, such as plant productivity, human disturbance, and the regional climate45,49. Our research indicated that the relationships of species diversity conformed to unimodal change patterns along various elevations for mountainous herb, shrub, and tree communities in a semi-arid region of the temperate zone. It can be concluded that vertical patterns of species diversity with a unimodal type may be a more universal phenomenon, relative to monotonic decreasing or increasing patterns of species diversity with increasing elevations.Factors influencing MDE at the plant community levelIn this study, more forbs and grasses were found at higher elevations, whereas more sedges occurred at lower elevations. The responses of the importance values of tree species to the altitudinal gradient demonstrated the following variation patterns: evergreen coniferous trees had higher importance values than deciduous coniferous trees, followed by deciduous broad-leaved trees. This showed that MDE was influenced by species life form. The species diversity of different life forms responded differently to the environment, and plant species with different life forms presented different diversity patterns along altitudinal gradients50. In New Zealand, the number of mountainous plants species decreased with increasing elevations and the total species number of all plants also decreased significantly, while species diversity had no significant distribution trend in response to elevation when plants with different life forms were considered under different layers of plant communities51.MDE is a common pattern of species diversity of mountainous plants with changing elevations. Our study area, located on the Loess Plateau of China, belongs to a semi-arid mountainous region in the temperate zone where the maximum species diversity of the tree community occurred at intermediate elevations. This finding was in accordance with the MDE hypothesis. Research from the Kinabalu Mountains in Sabah, Malaysia, indicated an obvious MDE pattern of species number linked to elevations52. On the Haleakala Mountains, Hawaii, USA, the highest species diversity occurred at intermediate elevations, which was also in accordance with the MDE hypothesis53. The MDE hypothesis was also proved by studies conducted in the Yu Mountains of Taiwan and the Emei Mountains of Szechwan in China54. However, the MDE pattern of species diversity in tree communities is caused by precipitation, which is the highest at intermediate elevations52. This situation also occurred in the herb community.There are many factors that affect the distribution of herbaceous plants, so the variations in species diversity in the herb community with elevations are complex55. In this study, we found that the herb community exhibited higher species diversity at higher elevations; more forbs and grasses were distributed at higher elevations, whereas more sedges were distributed at lower elevations. In the Siskiyou Mountains in Oregon, USA, the species diversity of herbaceous plants had a significantly positive correlation with elevation. This correlation occurred mainly due to an increase in the number of grass species, which was the primary reason that radiation was enhanced by a drastic reduction in community coverage as a result of increasing elevations. Consequently, there was an increase in the species diversity of herbaceous plants24. A decrease in species diversity with increasing elevations is a more familiar pattern for herbaceous plants in temperate56 and tropical22 forests.We also discovered that the family numbers of herb and shrub species all showed unimodal change patterns with high values at their central elevations in this semi-arid region. This conformed to the MDE hypothesis as well. In arid temperate grasslands, species diversity indicated an MDE distribution pattern. For example, the species diversity of herbaceous plants presented an MDE pattern in drought areas of the Siskiyou Mountains57. However, in semi-humid mountains in the temperate zone, the species diversity of herbaceous plants was principally in control of the community structure, and community coverage did not respond uniformly to elevation. Studies in New Zealand showed that there were no evident distribution trends for species diversity of herbaceous plants along elevations51. In low bush communities of Chile, the species diversity of herbaceous plants declined with increasing elevations after longstanding succession, but it increased during the early stage of succession58. Therefore, relationships between the species diversity of herb species and elevations were not completely clear in semi-humid regions.The major factors that control the distribution areas of species differ among different families and genera, and thus the vertical distribution patterns of species diversity differ22. We concluded that the family number of the tree species had a maximum distribution at lower elevations, unlike herb and shrub species; meanwhile, the responses of the importance values to the altitudinal gradient in the tree community were also different among evergreen coniferous trees, deciduous coniferous trees, and deciduous broad-leaved trees. These differences may have been related to environmental factors. Due to various distribution patterns of environmental parameters with elevations, the distribution patterns of species diversity showed large changes along elevations59. For example, the distribution of fern and Melastomataceae species is principally related to humidity, that of Acanthaceae and Bromeliaceae species is correlated with temperature, and that of Araceae species is related to transpiration59. Research conducted in the Gongga Mountains, China, showed that the species diversity with different floral components exhibited different distribution patterns along elevations due to differences in the environment and species origin60. We also discovered that the importance values of dominant families in the shrub (Rosaceae) and tree (Pinaceae) communities exhibited changing patterns in contrast to MDE. In our study, only the family numbers in the herb and shrub communities, as well as species diversity in the tree community, conformed to the MDE hypothesis. Therefore, we concluded that the MDE hypothesis of species diversity of mountainous vegetation is influenced by the species life form and family of different plant communities in the temperate semi-arid region of China.Factors influencing plant species diversity at the environment levelThe altitudinal distribution patterns of the plant community diversity had greater discrepancies in mountainous regions and between different community types, which might be connected to regional environmental conditions, relative heights of mountains, and the geological landscape35. Concerning the altitudes of mountains, serious human disturbances (e.g., deforestation, grazing, and land-use conversion) had negative effects on biodiversity in low-altitude regions61. In high-altitude regions, the cold climate slowed down plant growth and soil development, while other harsh environments exceeded the tolerance limits for growth of the majority of species, such as by intense solar radiation or large temperature differences between day and night62. In the middle-altitude regions, the species diversity was relatively higher due to less human disturbances and the formation of transition zones of plant species differentiation between the low- and high-altitude regions62. Hence, the plant community diversity and its altitudinal gradient patterns in mountainous regions were largely influenced by regional climate and human disturbances.Comparisons of the diversity at different levels indicated that the responses of the plant community diversity to the environment were not the same for diverse gradations, and different species exhibited different gradient patterns owing to restrictions from environmental factors63,64. The primary factors leading to the altitudinal differentiation of diversity included the temperature, moisture, soil nutrients, and succession process65. In our study area, compared with Guancen Mountain and Guandi Mountain, Wulu Mountain at the lower latitude of the Lvliang Mountains had a lower altitude and was located in the continental monsoon subhumid climate region of the warm temperate zone, making it suitable for the growth of secondary forests and shrub vegetation. However, the vegetation growth in the herb layer was restricted in Wulu Mountain, making diversity in the herb layer the greatest on Guandi Mountain at the middle latitude of the Lvliang Mountains35. This showed that the plant species diversity in the east of the Loess Plateau changed with the altitude, while being affected by complicated habitat conditions such as latitude and human disturbances. This characterization of the study area was correlated with the unimodal patterns observed. In this study, we observed that the family numbers of the herb and shrub communities presented unimodal patterns across an altitudinal gradient; the importance values of dominant families also presented unimodal patterns in the shrub and tree communities; and the species diversity indices of the herb, shrub, and tree communities conformed to unimodal change patterns following an altitudinal gradient as well.In our recent studies on the species diversity of herbaceous communities in the Lvliang Mountains66, we found that the results calculated for β-diversity using different indices revealed the highest value for the Cody index and the lowest value for the Bray–Curtis index at altitudes between 1900 and 2000 m, indicating that areas located between 1900 and 2000 m form a transition zone in which the herbaceous community undergoes a rapid process of species renewal and changes in its composition. The results for γ-diversity indicated a pattern of unimodal variation in relation to altitude. Changes in altitude gradient had highly significant impacts on changes in temperature and humidity, indicating that various environmental factors (notably humidity and temperature) and human disturbances had combined effects on changes in the values of the α-diversity indices.At present, it is widely believed that the formation of herbs in different life forms was principally impacted by precipitation, whereas in areas with similar rainfall, water, heat, and light conditions need to be considered. These conditions chiefly included average annual precipitation, accumulated temperature, and illumination time67. In this research, we observed that herb and tree species in different life forms showed different trends with altitudinal gradients in the Lvliang Mountains. At higher elevations, forbs and grasses grew well, whereas sedges grew well at lower elevations. The responses of different tree life forms to the altitudinal gradient were greater for evergreen coniferous tree species than for deciduous coniferous tree species and deciduous broad-leaved tree species. From the north to the south in the Lvliang Mountains, increases in the average annual precipitation increased the number of species and components of the annual herbs, while the hydrothermal matching requirements of Guandi Mountain at the middle latitude were preferred for annual herb growth in comparison with Guancen Mountain at the higher latitude67. However, considering whole mountains, the Lvliang Mountains located in the continental monsoon climate region of the warm temperate zone had four distinctive seasons with drought and wind in the spring, and a quick rise in air temperatures, and larger diurnal temperature differences35. These conditions conformed to the habitat features of herbs and trees. Hence, the hydrothermal distribution status affected by latitude and human disturbances determined the altitudinal distribution patterns of plant community diversity in the Lvliang Mountains.The MDE at different elevationsIn this study, we discussed the MDE of mountainous vegetation on the Loess Plateau with an elevational range from 1324 to 2745 m, including tree, shrub and herb community. This range was a very large elevation range for a case study, but a very short range in comparison to global elevation ranges, which extended from the sea level to well over 8000 m (though the highest locations did not have any vegetation). Therefore, owing to this limitation, the results we obtained in this research were suitable for lower elevation mountains in semi-arid areas.The MDE was changed with different elevations and vegetation layers. In studies on the Daiyun Mountains with an elevation from 900 to 1600 m, the phylogenetic diversity and species diversity of tree community indicated an intermediate high expansion pattern along elevations and their peak values all appeared at the elevation of 1200 m68. This conclusion conformed to the MDE pattern. In our studies on the Lvliang Mountains with an elevation from 1459 to 2610 m, higher species diversity of tree community was observed at intermediate elevations with a peak value being at the 2000 m, which conformed to the MDE pattern either. Therefore, at smaller elevations less than 2600 m, species diversity of tree conformed to the MDE pattern.When an elevation reached 2700 m on the Lvliang Mountains, the vegetation types changed to shrub and grass, and only their family numbers followed the MDE pattern across an altitudinal gradient. Slimily, the species richness of shrub and herb layer showed an obvious “lateral pattern” on an elevational gradient from 2950 to 3750 m on the Three River Headwater, both reaching the maximum value at the 3150 m; while with the rise of altitude, α diversity of shrub layer and herb layer showed a “wave”-shaped changed trend, reaching the lowest value at the 3550 m69. It illustrated that species diversity of shrub and grass did not conform to the MDE pattern completely at medium elevations from 2700 to 3700 m.As for an elevation extending from 3000 to 4400 m on the Gongga Mountains, the vegetation type was alpine meadow, and the species richness index presented an obvious unimodal pattern with a peak value at the 3850 m, which accorded with the MDE pattern70. Similarly, studies on an alpine meadow on the Gannan revealed that the number of richness, Shannon-Weiner index and phylogenetic diversity of plant community all showed a “humped-back” relationship with the increase of altitude from 3000 to 4000 m, and reached the maximum value at the 3800 m71,72. Thereby, at greater elevations more than 3800 m, species diversity of alpine meadow conformed to the MDE pattern.However, when an elevation exceeded 5000 m, research object on species diversity were not vegetation but animals along an altitudinal gradient. For example, on the Himalaya Mountains with an elevation from 3755 to 5016 m, the ant species richness illustrated a “unimodal curve” pattern along the rise of altitude, and the Shannon–Wiener index and Fisher α index of ant community commonly expressed the “Multi-Domain Effect” phenomenon73. Another research on mammalian richness was also conducted on the Himalaya Mountains. It concluded that most of elevational species richness patterns were hump-shaped from 100 to 6000 m on the Himalayas Mountains74. As a result, the MDE pattern was also extremely common in animal communities. More