Philippa Kaur

Tracing the symptomology of predation through macroscopic plaquesA culture bioassay (Expanded Microcoleus Mortality Assay, or EMMA) (Fig. 1 and see Materials and Methods) based on the capacity of a soil to induce complete mortality in the foundational biocrust cyanobacterium Microcoleus vaginatus helped us trace the pathogen detected in biocrust production facilities to the development of cm-sized plaques, or zones of cyanobacterial clearing, in natural biocrusts. These plaques were revealed to the naked eye (Fig. 2) when the soil was wet (i.e., after a rain event), as impacted areas would fail to green up by the migration of cyanobacteria to the surface21, enabling us to detect and quantify them with relative ease. Soil samples obtained from such plaques (n = 30) from different sites (n = 6; Table S1) in the US Southwest were invariably EMMA + , and the pathogens always filterable with pore sizes 0.45–1 µm but not larger, and always insensitive to the eukaryotic inhibitor cycloheximide, indicating the agent’s prokaryotic nature and small size, while paired samples from asymptomatic areas just outside the plaques were always EMMA- (Table S2). These end-point EMMA solutions never gave rise to cyanobacterial re-growth upon further incubation and maintained its infectivity of fresh cyanobacterial cultures for up to 6 months. A one-time, small-scale sampling across a plaque at intervals of 2 mm using microcoring22 showed that the boundary of the visible plaque demarcated exactly the end of infectivity, samples 0–2 mm outside the plaque proving non-infective. Further, inoculation of healthy, natural biocrusts with EMMA + suspensions resulted in the local development of biocrust plaques, and soil from these plaques was itself EMMA + , in partial fulfillment of Koch’s postulates. Yet, standard microbiological plating failed to yield any isolates that were EMMA + (we tested 30 unique isolates), even though standard plating with similar isolation efforts can successfully cultivate a large portion of heterotrophs from biocrusts23.Fig. 1: EMMA bioassay (Expanded Microcoleus Mortality Assay), used to study biocrust pathogens.a Typical visual progression of a positive EMMA inoculated with soil or culture to be tested, as used to test for pathogenicity to Microcoleus vaginatus PPC 9802 in the field and in enrichments. b Typical degradation of cyanobacterial biomass during an EMMA displayed through electron microscopy: healthy Microcoleus vaginatus PPC 9802 filaments (top) display abundant photosynthetic membranes (white arrows), peptidoglycan cross-walls (yellow arrows) and carboxysomes (green arrow). As infection proceeds (downwards), patent degradation of intracellular structures follows, leaving only cellular ghosts in the form of peptidoglycan wall remnants (yellow arrows), including the characteristically enlarged peptidoglycan “bumper” of terminal cells (red arrow). Intracellular bacilloid bacteria can sometimes be observed (blue arrow). Cyanobacterial cultures lose all viability. Scale bars = 1 µm. n = 250 images from 4 independent experiments. c Assay modification used in flow cytometry/cell sorting, showing enrichments positive for predation in the top two rows and those negative for predation below. d Test and controls in EMMA to ensure prokaryotic nature of the disease agent.Full size imageFig. 2: Symptomology in nature: biocrust plaques.Main: Macroscopic view of a soil surface colonized by cyanobacterial biocrusts and impacted by multiple plaques as taken after a rain in a quadrat used for field surveys. Insert: Close-up of a single plaque, showing well-demarcated boundaries and a typical central area of new cyanobacterial colonization.Full size imageCultivation, identification, and salient genomic traits of the cyanobacterial pathogenTo study these organisms, we turned to enrichment of pathogen/prey co-cultures based on repeated passages through EMMA and differential size filtration combined with dilution-to-extinction approaches, followed by purification with flow cytometry/cell sorting. The process was monitored by 16S rRNA gene amplicon sequencing, and eventually yielded a highly enriched co-culture of the cyanobacterium with a genetically homogenous (one single Amplicon Sequence Variant) population that made up more than 80% of reads (Fig. 3 a, b). We name the organism represented by this ASV Candidatus Cyanoraptor togatus. That it corresponds indeed to the predator is supported by the fact that of the 17 ASV’s detected in the final enrichment, only 10 were consistently detected at all infectious stages in the process and, among these, only our candidate ASV steadily increased in relative abundance through the enrichment process (Fig. 3 a, b). This final enrichment of C. togatus, LGM-1, constitutes the basis for downstream biological and molecular analyses. Its ASV was most similar to little-known members of the family Chitinophagaceae in the phylum Bacteroidetes. LGM-1’s genome was sequenced and assembled into a single 3.3 Mb contig with 1,781 putative and 1,328 hypothetical genes (Table S3), though most proteins had low identity (Fig. 4: Compiled paired ratios of functional parameters and compositional (relative) abundance in biocrusts across plaque boundaries (circles), red bars denoting the medians for each group of ratios, and bar background color denoting the p-values that the median is significantly different from unity (Wilcoxon paired ratio two-sided tests), where gray is non-significant (p  >  0.1), light orange is 0.05   > p   p  More

Hwang, H. Y. et al. Effect of cover cropping on the net global warming potential of rice paddy soil. Geoderma 292, 49–58 (2017).ADS 
CAS 
Article 

Google Scholar 
IPCC. Climate change 2013: The physical science basis. The Working Group I contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, 2013).
Google Scholar 
Cardoso, R. M., Soares, P. M. M., Lima, D. C. A. & Miranda, P. M. A. Mean and extreme temperatures in warming climate: EURO CORDEX and WRF regional climate high-resolution projection for Portugal. Clim. Dyn. 52, 129–157 (2019).Article 

Google Scholar 
Ding, T., Gao, H. & Li, W. J. Extreme high-temperature event in southern China in 2016 and the possible role of cross-equatorial flows. Int. J. Climatol. 38, 3579–3594 (2018).Article 

Google Scholar 
Escalas, A. et al. Functional diversity and redundancy across fish gut, sediment, and water bacterial communities. Environ. Microbiol. 19, 3268–3282 (2017).Article 

Google Scholar 
Philippot, L. et al. Loss in microbial diversity affects nitrogen cycling in soil. ISME J. 7, 1609–1619 (2013).CAS 
Article 

Google Scholar 
Li, Y. B. et al. Serratia spp. Are responsible for nitrogen fixation fueled by As(III) oxidation, a novel biogeochemical process identified in mine tailings. Environ. Sci. Technol 56, 2033–2043 (2022).ADS 
Article 

Google Scholar 
Jia, M., Gao, Z. W., Gu, H. J., Zhao, C. Y. & Han, G. D. Effects of precipitation change and nitrogen addition on the composition, diversity, and molecular ecological network of soil bacterial communities in a desert steppe. PLoS ONE 16, e0248194. https://doi.org/10.1371/journal.pone.0248194 (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Waghmode, T. R. et al. Response of nitrifier and denitrifier abundance and microbial community structure to experimental warming in an agricultural ecosystem. Front. Microbiol. 9, 474. https://doi.org/10.3389/fmicb.2018.00474 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Hu, Y. L., Wang, S., Niu, B., Chen, Q. & Zhang, G. Effect of increasing precipitation and warming on microbial community in Tibetan alpine steppe. Environ. Res. 189, 109917. https://doi.org/10.1016/j.envres.2020.109917 (2020).CAS 
Article 
PubMed 

Google Scholar 
Li, H. et al. Responses of soil bacterial communities to nitrogen deposition and precipitation increment are closely linked with aboveground community variation. Microb. Ecol. 71, 974–989 (2016).CAS 
Article 

Google Scholar 
Wang, H. et al. Experimental warming reduced topsoil carbon content and increased soil bacterial diversity in a subtropical planted forest. Soil Biol. Biochem. 133, 155–164 (2019).CAS 
Article 

Google Scholar 
Haumann, F. A., Gruber, N. & Münnich, M. Sea-Ice Induced Southern Ocean Subsurface Warming and Surface Cooling in a Warming Climate. AGU Advances 1, e2019AV000132. https://doi.org/10.1029/2019AV000132 (2020).ADS 
Article 

Google Scholar 
Ji, F., Wu, Z. H., Huang, J. P. & Chassignet, E. P. Evolution of land surface air temperature trend. Nat. Clim. Chang. 4, 462–466 (2014).ADS 
Article 

Google Scholar 
Sabri, N. S. A., Zakaria, Z., Mohamad, S. E., Jaafar, A. B. & Hara, H. Importance of soil temperature for the growth of temperate crops under a tropical climate and functional role of soil microbial diversity. Microbes Environ. 33, 144–150 (2018).Article 

Google Scholar 
McGrady-Steed, J. & Morin, P. T. Biodiversity, density compensation, and the dynamics of populations and functional groups. Ecology 81, 361–373 (2000).Article 

Google Scholar 
Jiang, L. Density compensation can cause no effect of biodiversity on ecosystem function. Oikos 116, 324–334 (2007).Article 

Google Scholar 
Faust, K. & Raes, J. Microbial interactions: From networks to models. Nat. Rev. Microbiol. 10, 538. https://doi.org/10.1038/nrmicro283 (2012).CAS 
Article 
PubMed 

Google Scholar 
Gao, X. X. et al. Revegetation significantly increased the bacterial-fungal interactions in different successional stages of alpine grasslands on the Qinghai-Tibetan Plateau. CATENA 205, 105385. https://doi.org/10.1016/j.catena.2021.105385 (2021).CAS 
Article 

Google Scholar 
Morriën, E. et al. Soil networks become more connected and take up more carbon as nature restoration progresses. Nat. Commun. 8, 14349. https://doi.org/10.1038/ncomms14349 (2017).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Banerjee, S. et al. Agricultural intensification reduces microbial network complexity and the abundance of keystone taxa in roots. ISME J. 13, 1722–1736 (2019).Article 

Google Scholar 
Pržulj, N. & Malod-Dognin, N. Network analytics in the age of big data. Science 353, 123–124 (2016).ADS 
Article 

Google Scholar 
Ratzke, C., Barrere, J. M. R. & Gore, J. Strength of species interactions determines biodiversity and stability in microbial communities. Nat. Ecol. Evol. 4, 376–383 (2020).Article 

Google Scholar 
Fuhrman, J. A. Microbial community structure and its functional implications. Nature 45, 193–199 (2009).ADS 
Article 

Google Scholar 
Zhao, M. X., Cong, J., Cheng, J. M., Qi, Q. & Zhang, Y. G. Soil microbial community assembly and interactions are constrained by nitrogen and phosphorus in broadleaf forests of southern China. Forest 11, 285. https://doi.org/10.3390/f11030285 (2020).Article 

Google Scholar 
Wan, X. L. et al. Biogeographic patterns of microbial association networks in paddy soil within Eastern China. Soil Biol. Biochem. 142, 07696. https://doi.org/10.1016/j.soilbio.2019.107696 (2020).CAS 
Article 

Google Scholar 
Yuan, M. M., Guo, X., Wu, L., Zhang, Y. & Zhou, J. Climate warming enhances microbial network complexity and stability. Nat. Clim. Change 11, 343–348 (2021).ADS 
Article 

Google Scholar 
Lassaletta, L. et al. Food and feed trade as a driver in the global nitrogen cycle: 50-year trends. Biogeochemistry 11, 225–241 (2014).Article 

Google Scholar 
Phoenix, G. K. et al. Impacts of atmospheric nitrogen deposition: Responses of multiple plant and soil parameters across contrasting ecosystems in long-term field experiments. Glob. Change Biol. 18, 1197–1215 (2012).ADS 
Article 

Google Scholar 
Nakaji, T., Fukami, M., Dokiya, Y. & Izuta, T. Effects of high nitrogen load on growth, photosynthesis and nutrient status of Cryptomeria japonica and Pinus densiflora seedlings. Trees-Struct. Funct. 15, 453–461 (2001).CAS 
Article 

Google Scholar 
Wang, H. Y. et al. Reduction in nitrogen fertilizer use results in increased rice yields and improved environmental protection. Int. J. Agric. Sustain. 15, 681–692 (2017).Article 

Google Scholar 
Zhou, X. G. & Wu, F. Z. Land-use conversion from open field to greenhouse cultivation differently affected the diversities and assembly processes of soil abundant and rare fungal communities. Sci. Total Environ. 788, 147751. https://doi.org/10.1016/j.scitotenv.2021.147751 (2021).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Guo, H. et al. Long-term nitrogen & phosphorus additions reduce soil microbial respiration but increase its temperature sensitivity in a Tibetan alpine meadow. Soil Biol. Biochem. 113, 26–34 (2017).CAS 
Article 

Google Scholar 
Zhang, C. et al. Effects of simulated nitrogen deposition on soil respiration components and their temperature sensitivities in a semiarid grassland. Soil Biol. Biochem. 75, 113–123 (2014).CAS 
Article 

Google Scholar 
Zhang, J. J. et al. Different responses of soil respiration and its components to nitrogen and phosphorus addition in a subtropical secondary forest. For. Ecosyst. 8, 37. https://doi.org/10.1186/s40663-021-00313-z (2021).Article 

Google Scholar 
Norse, D. & Ju, X. T. Environmental costs of China’s food security. Agric. Ecosyst. Environ. 209, 5–14 (2015).Article 

Google Scholar 
Xu, H. F., Du, H., Zeng, F. P., Song, T. Q. & Peng, W. X. Diminished rhizosphere and bulk soil microbial abundance and diversity across succession stages in Karst area, southwest China. Appl. Soil Ecol. 158, 103799. https://doi.org/10.1016/j.apsoil.2020.103799 (2020).Article 

Google Scholar 
Li, Y. B. et al. Arsenic and antimony co-contamination influences on soil microbial community composition and functions: Relevance to arsenic resistance and carbon, nitrogen, and sulfur cycling. Environ. Int. 153, 106522. https://doi.org/10.1016/j.envint.2021.106522 (2021).CAS 
Article 
PubMed 

Google Scholar 
Zhou, J. & Fong, J. J. Strong agricultural management effects on soil microbial community in a non-experimental agroecosystem. Appl. Soil Ecol. 165, 103970. https://doi.org/10.1016/j.apsoil.2021.103970 (2021).Article 

Google Scholar 
Bárcenas-Moreno, G., Gómez-Brandón, M., Rousk, J. & Bååth, E. Adaptation of soil microbial communities to temperature: Comparison of fungi and bacteria in a laboratory experiment. Glob. Chang. Biol. 15, 2950–2957 (2009).ADS 
Article 

Google Scholar 
Tan, E. H., Zou, W., Zheng, Z., Yan, X. & Kao, S. J. Warming stimulates sediment denitrification at the expense of anaerobic ammonium oxidation. Nat. Clim. Change 10, 349–355 (2020).ADS 
CAS 
Article 

Google Scholar 
Supramaniam, Y., Chong, C. W., Silvaraj, S. & Tan, K. P. Effect of short term variation in temperature and water content on the bacterial community in a tropical soil. Appl Soil Ecol. 107, 279–289 (2016).Article 

Google Scholar 
Zhu, Y. Z., Li, Y. Y., Zheng, N. G., Chapman, S. J. & Yao, H. Y. Similar but not identical resuscitation trajectories of the soil microbial community based on either DNA or RNA after flooding. Agronomy 10, 502. https://doi.org/10.3390/agronomy10040502 (2020).CAS 
Article 

Google Scholar 
Donhauser, J., Qi, W., Bergk-Pinto, B. & Frey, B. High temperatures enhance the microbial genetic potential to recycle C and N from necromass in high-mountain soils. Glob. Chang. Biol. 27, 1365–1386 (2021).ADS 
Article 

Google Scholar 
Santoyo, G., Hernandez-Pacheco, C., Hernandez-Salmeron, J. & Hernandez-Leon, R. The role of abiotic factors modulating the plant-microbe-soil interactions: Toward sustainable agriculture. A review. Span. J. Agric. Res. 15, e03R01-e11. https://doi.org/10.5424/sjar/2017151-9990 (2017).Article 

Google Scholar 
Lefcheck, J. S. et al. Biodiversity enhances ecosystem multifunctionality across trophic levels and habitats. Nat. Commun. 6, 6936. https://doi.org/10.1038/ncomms7936 (2015).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Cardinale, B. J. et al. Corrigendum: Biodiversity loss and its impact on humanity. Nature 486, 59–67 (2012).ADS 
CAS 
Article 

Google Scholar 
Ma, B., Wang, H., Dsouza, M., Lou, J. & Xu, J. Geographic patterns of co-occurrence network topological features for soil microbiota at continental scale in eastern China. ISME J. 10, 1891–1901 (2016).CAS 
Article 

Google Scholar 
Trivedi, C. et al. Losses in microbial functional diversity reduce the rate of key soil processes. Soil Biol. Biochem. 135, 267–274 (2019).CAS 
Article 

Google Scholar 
Melanie, K. et al. Effects of season and experimental warming on the bacterial community in a temperate mountain forest soil assessed by 16S rRNA gene pyrosequencing. FEMS Microbiol. Ecol. 82, 551–562 (2012).Article 

Google Scholar 
Zheng, H. F., Liu, Y., Chen, Y., Zhang, J. & Chen, Q. Short-term warming shifts microbial nutrient limitation without changing the bacterial community structure in an alpine timberline of the eastern Tibetan Plateau. Geoderma 360, 113985. https://doi.org/10.1016/j.geoderma.2019.113985 (2020).ADS 
CAS 
Article 

Google Scholar 
Finlay, B. J. & Cooper, J. L. Microbial diversity and ecosystem function. CEH Integrating Fund second progress report to the Director, Centre for Ecology and Hydrology Nov 1996–Sept (1997).Xing, X. Y. et al. Warming shapes nirS- and nosZ-type denitrifier communities and stimulates N2O emission in acidic paddy soil. Appl. Environ. Microbiol. 87, e02965-e3020. https://doi.org/10.1128/AEM.0296520 (2021).CAS 
Article 
PubMed Central 

Google Scholar 
Lin, Y. T., Whitman, W. B., Coleman, D. C., Jien, S. H. & Chiu, C. Y. Soil bacterial communities at the treeline in subtropical alpine areas. CATENA 201, 105205. https://doi.org/10.1016/j.catena.2021.105205 (2021).CAS 
Article 

Google Scholar 
Wang, J. C. et al. Impacts of inorganic and organic fertilization treatments on bacterial and fungal communities in a paddy soil. Appl. Soil Ecol. 112, 42–50 (2017).Article 

Google Scholar 
Chacón, J. M., Shaw, A. K. & Harcombe, W. R. Increasing growth rate slows adaptation when genotypes compete for diffusing resources. PLoS Comput. Biol. 16, e1007585. https://doi.org/10.1371/journal.pcbi.1007585 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Hartley, I. P., Hopkins, D. W., Garnett, M. H., Sommerkorn, M. & Wookey, P. A. Soil microbial respiration in arctic soil does not acclimate to temperature. Ecol. Lett. 11, 1092–1100 (2008).Article 

Google Scholar 
Baath, E. Growth rates of bacterial communities in soils at varying pH: A comparison of the thymidine and leucine incorporation techniques. Microb. Ecol. 36, 316–327 (1998).CAS 
Article 

Google Scholar 
Qin, H. L. et al. Soil moisture and activity of nitrite- and nitrous oxide-reducing microbes enhanced nitrous oxide emissions in fallow paddy soils. Biol. Fertil. Soils 56, 53–67 (2020).CAS 
Article 

Google Scholar 
Chen, Z. et al. Impact of long term fertilization on the composition of denitrifier communities based on nitrite reductase analyses in a paddy soil. Microb. Ecol. 60, 850–861 (2010).CAS 
Article 

Google Scholar 
Wei, G. S. et al. Similar drivers but different effects lead to distinct ecological patterns of soil bacterial and archaeal communities. Soil Biol. Biochem. 144, 107759. https://doi.org/10.1016/j.soilbio.2020.107759 (2020).CAS 
Article 

Google Scholar 
Bastian, F., Bouziri, L., Nicolardot, B. & Ranjard, A. L. Impact of wheat straw decomposition on successional patterns of soil microbial community structure. Soil Biol. Biochem. 41, 262–275 (2009).CAS 
Article 

Google Scholar 
Levins, R. Evolution in Changing Environments: Some Theoretical Explorations (Princeton University Press, 1968).Book 

Google Scholar  More

Genetic structure of the parental populationBased on the lnPD and ΔK values obtained using STRUCTURE, we identified two genetic groups within the DHS Acer population (Supplementary Fig. S1). The q value from STRUCTURE analysis represents the proportion of ancestral origin28 (Fig. 2a). Among the 70 individual trees, 72.9% were assigned a q value smaller than 0.1 or larger than 0.9, thereby signifying a typical bimodal distribution (Fig. 2b). Individuals with q value greater than 0.9 and consistent genetic origin from the NEA region were defined as the NEA lineage (hereafter “NEA-DHS”), whereas those with values less than 0.1 and with consistent genetic origin from the SEA region were defined as the SEA lineage (hereafter “SEA-DHS”). Individuals with intermediate q value between 0.1 and 0.9 were defined as hybrid genetic types (hereafter “Hybrid-DHS”). Accordingly, we identified 27 SEA-DHS (38.6%), 24 NEA-DHS (34.3%), and 19 Hybrid-DHS (27.1%) (Fig. 2b).Figure 2Genetic structure of the parental and offspring population. (a) Bar plots illustrating the genetic composition of the adult (leaf) and offspring (fruit) populations in the Daheishan National Nature Reserve (DHS). Each individual is represented by a line partitioned into color segments corresponding to its ancestral proportion. Red color represents the ancestral proportion of Southern East Asia lineage. Green color represents the ancestral proportion of Northern East Asia lineage. Black lines in bar plots of leaf population separate individuals with ancestral proportion (q value) bigger than 0.9 or smaller than 0.1 from hybrids (0.1  0.5) produced by the SEA-DHS were obtained from a single tree, which was identified as SEA-DHS based on the DHS-only dataset, although it was indicated to be Hybrid-DHS based on the whole-range dataset. The Hybrid-DHS maternal trees produced 17.6% pure SEA-DHS seeds, 57.6% pure NEA-DHS seeds, and 24.7% hybrid seeds.Flowering phenologyThe sexual system of Acer has four phenotypes: duodichogamous, protogynous, protandrous, and male31. Hence, there are three functional sex types: (1) “Male I” flowers open earlier than “Female” flowers, with mature stamens, no style, and ovary; (2) “Female” flowers have mature pistils, short filaments, and indehiscence anthers; (3) “Male II” flowers open later than “Female” flowers, with mature stamens, ovaries, and separated stigmas. Duodichogamy is characterized by “Male I,” “Female,” and “Male II” types; protandry by “Male I” and “Female” types; and protogyny by “Female” and “Male II” types31.During the flowering season, we monitored a total of 10,074 flowers produced by 29 trees (Fig. 2d), among which one tree (SEA-DHS) was protandrous, four trees (three Hybrid-DHS and one NEA-DHS) were protogynous, and the remaining 24 trees were duodichogamous. We observed that the blooming phenology of SEA-DHS and NEA-DHS differed significantly to most assessed phenological indices, with a single exception being a marginally significant difference in the peak blooming time of Male I (Table 1). Compared with NEA-DHS, SEA-DHS were characterized by significantly later flowering phenology, with Male I commencement and cessation of blooming being on average two and three days later, respectively. Similarly, the commencement, peak, and cessation of Female occurred later by averages of 4, 4, and 5 days, respectively, whereas those of Male II occurred later by 5, 4, and 5 days, respectively. Furthermore, the duration of blooming was significantly longer in the SEA-DHS group than in the NEA-DHS group by three days. In the case of Hybrid-DHS, the values obtained for all assessed phenological indices were intermediate between those of the two parental types. Among these, the values of the six indices differed significantly from one or the other parental types, with the majority (5/6) differing from those of the SEA-DHS. Thus, phenologically, Hybrid-DHS appeared to be closer to NEA-DHS.Table 1 Flowering phenology of SEA-DHS, Hybrid-DHS, and NEA-DHS.Full size tableHowever, despite the differing phenology of the SEA-DHS and NEA-DHS, we observed instances of overlap in the blooming periods of male or female flowers in one genetic type with those of flowers of the opposite sex in another genetic type. For example, the peak of Female among NEA-DHS (11.67 ± 0.67) was found to coincide with the peak of Male I (11.44 ± 1.06; p = 0.879) in SEA-DHS. Similarly, Female blooming in the SEA-DHS peaked (16.11 ± 1.09) just 1 d after the peak of Male II (15.50 ± 0.43) in the NEA-DHS (p = 0.667), which at this time still retained an abundance of male flowers in bloom. In contrast, we detected no overlapping phenology with respect to the blooming of Male I of NEA-DHS or Male II of SEA-DHS with the Female in another genetic type.Morphological variation of leaves and fruitLeaves Among the eight leaf indices, all except InfectionRatio were significantly different between lineages. Generally, the leaves of NEA-DHS were found to have seven lobes, whereas those of SEA-DHS were typically five lobed (Lobes#), thereby contributing to significantly larger leaves in NEA-DHS than in SEA-DHS (TotalArea). Furthermore, NEA-DHS leaves had shorter and wider central lobes (CentralLength and CentralWidth), as well as an earlier and narrower inflection of the central lobes (InflectionLength and InflectionWidth), compared with those of SEA-DHS (Table 2). Six indices had correlation coefficients of less than 0.7, which were used for principal component analysis (PCA) analysis (Supplementary Table S2). The first two axes of the PCA were found to explain 63.7% of the variation in leaf morphology (Fig. 3a), with InflectionLength, CentralLength, and CentralRatio contributing the most to the first axis (38.2%), whereas TotalArea contributed the most to the second axis (25.5%) (Supplementary Table S3). The leaves of SEA-DHS and NEA-DHS plants were largely clustered in separate groups (Fig. 3a). However, all indices were continuous variables with large overlaps between the lineages (Table 2). For example, NEA-DHS had a significantly larger leaf area (21.06–88.70 cm2) than SEA-DHS (11.34–70.09 cm2). The shape of the central lobe is another major leaf trait that distinguishes between the two species. NEA-DHS had a shorter and wider central lobe (CentralRatio:0.67–2.49), while SEA-DHS had a longer and narrower central lobe (CentralRatio:0.9–3.46).Table 2 Morphological variation in the leaves and fruits of Acer trees in the Daheishan National Nature Reserve.Full size tableFigure 3Morphological variation in the leaves (a) and fruits (b) of southern and northern East Asia lineages of the Acer species complex in the Daheishan National Nature Reserve based on principal component analysis. SEA-DHS: Southern East Asia lineage of the Acer species complex in the DHS; NEA-DHS: Northern East Asia lineage of the Acer species complex in the DHS; Hybrid-DHS: hybrids between SEA-DHS and NEA-DHS lineages.Full size imageWith regard to Hybrid-DHS, the leaves were morphologically intermediate between those of the two parental types (Fig. 3a), as were the values of the assessed morphological trait indices (Table 2).Fruits 11 indices of fruits were significantly different between lineages. NEA-DHS tend to be characterized by smaller fruits (FruitLength and FruitWidth), seeds (SeedLength, SeedWidth and JunctionWidth), and fruit wings (WingLength and WingWidth). Moreover, the seed wings of NEA-DHS fruits are typically oriented at an obtuse angle, whereas those of SEA-DHS fruits tend to be aligned at a right angle (FruitAngle). The length ratio of the wing and seed (Wing:Seed) was larger in NEA-DHS than in SEA-DHS (1.24 vs 1.06, respectively, Table 2). Eight indices had correlation coefficients of less than 0.7, which were retained for PCA analysis (Supplementary Table S4). The first two axes of the PCA explained 58.4% of the variation in fruit morphology (Fig. 3b), with JunctionWidth and SeedLength contributing the most to the first axis (35.1%), whereas SeedRatio and WingRatio contributed the most to the second axis (23.3%) (Supplementary Table S3). The fruits of SEA-DHS and NEA-DHS plants were largely clustered in separate groups, with most fruits of SEA-DHS having negative values in Axis 1, while most fruits of NEA-DHS having positive values (Fig. 3b). Both JunctionWidth and SeedLength in Axis 1 reflect the size of the seed. NEA-DHS had smaller seed (SeedLength: 0.63–1.21 cm, SeedWidth:0.43–0.75 cm), while larger seed in SEA-DHS (SeedLength:0.79–1.49 cm, SeedWidth:0.49–0.93 cm). All indices were continuous variables with large overlaps between the lineages (Table 2).The morphology of Hybrid-DHS fruits was generally intermediate between that of the two parental types (Fig. 3b), as reflected in the values of the different morphological traits. The exceptions in this regard were FruitLength, WingLength, as well as two ratio indices (SeedRatio and WingRatio), with hybrid trees typically producing longer fruit with longer fruit wings (Table 2).Ecological niche divergence between NEA and SEAWe found a positive correlation between q value from Structure analysis and altitude (Pearson’s r = 0.83, p  670 m), whereas SEA-DHS was clustered at the foothill ( More

IPBES. Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. Zenodo https://doi.org/10.5281/zenodo.5657041 (2019).Laurance, W. F. et al. A global strategy for road building. Nature 513, 229–232 (2014).ADS 
CAS 
Article 

Google Scholar 
Ibisch, P. L. et al. A global map of roadless areas and their conservation status. Science 354, 1423–1427 (2016).ADS 
CAS 
Article 

Google Scholar 
Forman, R. T. T., Sperling, D. & Bissonette, J. A. Road Ecology: Science and Solutions. (Island Pr, 2003).van der Ree, R., Smith, D. J. & Grilo, C. Handbook of Road Ecology. (John Wiley & Sons, 2015).Laender. Hunting Statistics. Game casualties 2017/2018: furred game (red deer, roe deer, chamois, moufflon) https://www.statistik.at/web_en/statistics/Economy/agriculture_and_forestry/livestock_animal_production/hunting/index.html (2018).Steiner, W., Leisch, F. & Hacklander, K. A review on the temporal pattern of deer-vehicle accidents: Impact of seasonal, diurnal and lunar effects in cervids. Accident; analysis and prevention 66, (2014).Kioko, J. et al. Driver knowledge and attitudes on animal vehicle collisions in Northern Tanzania. TROPICAL CONSERVATION SCIENCE 8, 352–366 (2015).Article 

Google Scholar 
Bíl, M., Andrášik, R. & Janoška, Z. Identification of hazardous road locations of traffic accidents by means of kernel density estimation and cluster significance evaluation. Accident Analysis & Prevention 55, 265–273 (2013).Article 

Google Scholar 
Page, Y. A statistical model to compare road mortality in OECD countries. Accident Analysis and Prevention 33, 371–385 (2001).CAS 
Article 

Google Scholar 
Teixeira, F. Z. et al. Are Road-kill Hotspots Coincident among Different Vertebrate Groups? Oecologia Australis 17, 36–47 (2017).Article 

Google Scholar 
Canova, L. & Balestrieri, A. Long-term monitoring by roadkill counts of mammal populations living in intensively cultivated landscapes. Biodivers Conserv https://doi.org/10.1007/s10531-018-1638-3 (2018).Brehme, C. S., Hathaway, S. A. & Fisher, R. N. An objective road risk assessment method for multiple species: ranking 166 reptiles and amphibians in California. Landscape Ecol 33, 911–935 (2018).Article 

Google Scholar 
Heigl, F. et al. Comparing Road-Kill Datasets from Hunters and Citizen Scientists in a Landscape Context. Remote Sensing 8, (2016).Heigl, F., Horvath, K., Laaha, G. & Zaller, J. G. Amphibian and reptile road-kills on tertiary roads in relation to landscape structure: using a citizen science approach with open-access land cover data. BMC Ecol 17, 24 (2017).Article 

Google Scholar 
Dörler, D. & Heigl, F. A decrease in reports on road-killed animals based on citizen science during COVID-19 lockdown. PeerJ 9, e12464 (2021).Article 

Google Scholar 
Peer, M. et al. Predicting spring migration of two European amphibian species with plant phenology using citizen science data. Sci Rep 11, 21611 (2021).ADS 
CAS 
Article 

Google Scholar 
Schwartz, A. L. W. UK Roadkill Records. The Global Biodiversity Information Facility https://doi.org/10.15468/r3xakd (2018).Lin, T. The Taiwan Roadkill Observation Network Data Set. Version 1.3. The Global Biodiversity Information Facility https://doi.org/10.15468/cidkqi (2018).Chandler, M. et al. Contribution of citizen science towards international biodiversity monitoring. Biological Conservation 213, 280–294 (2017).Article 

Google Scholar 
Périquet, S., Roxburgh, L., le Roux, A. & Collinson, W. J. Testing the Value of Citizen Science for Roadkill Studies: A Case Study from South Africa. Front. Ecol. Evol. 6, (2018).Abra, F. D., Huijser, M. P., Pereira, C. S. & Ferraz, K. M. P. M. B. How reliable are your data? Verifying species identification of road-killed mammals recorded by road maintenance personnel in São Paulo State, Brazil. Biological Conservation 225, 42–52 (2018).Article 

Google Scholar 
Bíl, M., Kubeček, J., Sedoník, J. & Andrášik, R. Srazenazver.cz: A system for evidence of animal-vehicle collisions along transportation networks. Biological Conservation 213, 167–174 (2017). Part A.Article 

Google Scholar 
Vercayie, D. & Herremans, M. Citizen science and smartphones take roadkill monitoring to the next level. Nature Conservation 11, 29–40 (2015).Article 

Google Scholar 
Waetjen, D. P. & Shilling, F. M. Large Extent Volunteer Roadkill and Wildlife Observation Systems as Sources of Reliable Data. Front. Ecol. Evol. 5, (2017).Shilling, F. M., Perkins, S. E. & Collinson, W. Wildlife/Roadkill Observation and Reporting Systems. in Handbook of Road Ecology 492–501 (John Wiley & Sons, 2015).Eitzel, M. V. et al. Citizen Science Terminology Matters: Exploring Key Terms. Citizen Science: Theory and Practice 2, 1–20 (2017).
Google Scholar 
Haklay, M. et al. Contours of citizen science: a vignette study. Royal Society Open Science 8, 202108 (2021).ADS 
Article 

Google Scholar 
Heigl, F., Kieslinger, B., Paul, K. T., Uhlik, J. & Dörler, D. Opinion: Toward an international definition of citizen science. PNAS 116, 8089–8092 (2019).CAS 
Article 

Google Scholar 
Heigl, F. et al. Quality Criteria for Citizen Science Projects on Österreich forscht | Version 1.1. Open Science Framework https://doi.org/10.17605/OSF.IO/48J27 (2018).Heigl, F. et al. Co-Creating and Implementing Quality Criteria for Citizen Science. Citizen Science: Theory and Practice 5, 23 (2020).
Google Scholar 
Heigl, F. & Zaller, J. G. Using a Citizen Science Approach in Higher Education: a Case Study reporting Roadkills in Austria. Human Computation 1, (2014).University of Natural Resources and Life Sciences, Vienna. Roadkill, The Global Biodiversity Information Facility, https://doi.org/10.15468/ejb47y (2021).Heigl, F. & Roadkill Community. Roadkill Dataset 2014-2020 Quality level 2, Zenodo, https://doi.org/10.5281/zenodo.5878813 (2022).August, T. A. et al. Citizen meets social science: predicting volunteer involvement in a global freshwater monitoring experiment. Freshwater Science 38, 321–331 (2019).Article 

Google Scholar 
IUCN. The IUCN Red List of Threatened Species. Version 2021-3. IUCN Red List of Threatened Species https://www.iucnredlist.org/en (2021). More

Shotgun sequencing and generation of metagenome-assembled genomesIn our previous study, 78 Holstein Friesian dairy cows were sampled for rumen content, metagenomic shotgun sequencing was carried out, and raw Illumina sequencing reads were assembled into contigs using megahit assembler using default settings [7]. We used a pooled assembly of the original 78 samples to increase the quality of the metagenome-assembled genomes (MAGs) with the syntax: megahit [14] -t 60 -m 0.5 −1 [Illumina R1 files] −2 [Illumina R2 files]. Next, the assembled contigs were indexed using BBMap [15]: bbmap.sh threads = 60 ref = [contigs filename]. Thereafter, reads from each sample were mapped to the assembled contigs using BBTools’ bbwrap.sh script. In order to determine the depth (coverage) of each contig within each sample, the gi_summarize_bam_contig_depths tool was applied with the parameters: gi_summarize_bam_contig_depths –outputDepth depth.txt –pairedContigs paired.txt *.bam –outputDepth depth.txt –pairedContigs paired.txt.Using the depth information, metabat2 [16] was executed to bind genes together into reconstructed genomes, with parameters: metabat2 -t40 -a depth.txt.To evaluate genomic bin quality, we used the CheckM [17] tool, with parameters: checkm lineage_wf [in directory] [out directory] -x faa –genes -t10.Preparing proteomic search libraryWe generated 93 unique high-quality MAGs, and further increased our MAG database by including phyla that were not represented in our set of MAGs. In order to do so, we used the published compendium of 4,941 rumen metagenome-assembled genomes [18] and dereplicated those MAGs using dRep [19]. We then selected MAGs from phylum Spirochaetes, Actinomycetota, Proteobacteria, Firmicutes, Elusimicrobia, Bacillota, Fibrobacteres and Fusobacteria, which had the highest mean coverage in our samples as calculated using BBMap and gi_summarize_bam_contig_depths as described above [15]. This strategy minimized the false discovery rate (FDR), that would have been obtained if larger and unspecific databases would have been employed [20] and allowed the addition of 14 MAGs to our database.In order to create the proteomic search library, genes were identified along the 107 MAGs using the Prodigal tool [21], with parameters: prodigal meta and translated in silico into proteins, using the same tool. Replicates sequences were removed. Protein sequences from the hosting animal (Bos taurus) and common contaminant protein sequences (64,701 in total) were added to the proteomic search library in order to avoid erroneous target protein identification originating from the host or common contaminants. Finally, in order to subsequently assess the percentage of false-positive identifications within the proteomic search [22], the proteomic search library sequences were reversed in order and served as a decoy database.Proteomic analysisThe bacterial fraction from rumen fluid of the 12 selected animals selected from extreme feed efficiency phenotypes, were obtained at the same time as the samples analyzed for metagenomics and stored at −20 °C until extraction. To extract total proteins, a modified protocol from Deusch and Seifert was used [23]. Briefly, cell pellets were resuspended in 100 µl in 50 mM Tris-HCl (pH 7.5; 0.1 mg/ml chloramphenicol; 1 mM phenylmethylsulfonyl fluoride (PMSF)) and incubated for 10 min at 60 °C and 1200 rpm in a thermo-mixer after addition of 150 µl 20 mM Tris-HCl (pH 7.5; 2% sodium dodecyl sulfate (SDS)). After the addition of 500 µl DNAse buffer (20 mM Tris-HCl pH 7.5; 0.1 mg/ml MgCl2, 1 mM PMSF, 1 μg/ml DNAse I), the cells were lysed by ultra-sonication (amplitude 51–60%; cycle 0.5; 4 × 2 min) on ice, incubated in the thermo-mixer (10 min at 37 °C and 1,200 rpm) and centrifuged at 10,000 × g for 10 min at 4 °C. The supernatant was collected and centrifuged again. The proteins in the supernatant were precipitated by adding 20% pre-cooled trichloroacetic acid (TCA; 20% v/v). After centrifugation (12,000 × g; 30 min; 4 °C), the protein pellets were washed twice in pre-cooled (−20 °C) acetone (2 × 10 min; 12,000 × g; 4 °C) and dried by vacuum centrifugation. The protein pellet was resuspended in 2× SDS sample buffer (4% SDS (w/v); 20% glycerin (w/v); 100 mM Tris-HCl pH 6.8; a pinch of bromophenol blue, 3.6% 2‑mercaptoethanol (v/v)) by 5 min sonication bath and vortexing. Samples were incubated for 5 min at 95 °C and separated by 1D SDS-PAGE (Criterion TG 4-20% Precast Midi Gel, BIO-RAD Laboratories, Inc., USA).As previously described, after fixation and staining, each gel line was cut into 10 pieces, destained, desiccated, and rehydrated in trypsin [24]. The in-gel digest was performed by incubation overnight at 37 °C. Peptides were eluted with Aq. dest. by sonication for 15 min The sample volume was reduced in a vacuum centrifuge.Before MS analysis, the tryptic peptide mixture was loaded on an Easy-nLC II or Easy-nLC 1000 (Thermo Fisher Scientific, USA) system equipped with an in-house built 20 cm column (inner diameter 100 µm; outer diameter 360 µm) filled with ReproSil-Pur 120 C18-AQ reversed-phase material (3 µm particles, Dr. Maisch GmbH, Germany). Peptides were eluted with a nonlinear 156 min gradient from 1 to 99% solvent B (95% acetonitrile (v/v); 0.1% acetic acid (v/v)) in solvent A (0.1% acetic acid (v/v)) with a flow rate of 300 ml/min and injected online into an LTQ Orbitrap Velos or Orbitrap Velos Pro (Thermo Fisher Scientific, USA). Overview scan at a resolution of 30,000 in the Orbitrap in a range of 300-2,000 m/z was followed by 20 MS/MS fragment scans of the 20 most abundant precursor ions. Ions without detected charge state as well as singly charged ions were excluded from MS/MS analysis. Original raw spectra files were converted into the common mzXML format, in order to further process it in downstream analysis. The spectra file from each proteomic run of a given sample was searched against the protein search library, using the Comet [25] search engine with default settings.The TPP pipeline (Trans Proteomic Pipeline) [26] was used to further process the Comet [25, 27] search results and produce a protein abundance table for each sample. In detail, PeptideProphet [28] was applied to validate peptide assignments, with filtering criteria set to probability of 0.001, accurate mass binning, non-parametric errors model (decoy model) and decoy hits reporting. In addition, iProphet [28, 29] was applied to refine peptide identifications coming from PeptideProphet. Finally, ProteinProphet [28,29,30] was applied to statistically validate peptide identifications at the protein level. This was carried out using the command: xinteract -N[my_sample_nick].pep.xml -THREADS = 40 -p0.001 -l6 -PPM -OAPd -dREVERSE_ -ip [file1].pep.xml [file2].pep.xml.. [fileN].pep.xml  > xinteract.out 2  > xinteract.err. Then, TPP GUI was used in order to produce a protein table from the resulting ProtXML files (extension ipro.prot.xml).Subsequently, proteins that had an identification probability < 0.9 were also removed as well as proteins supported with less than 2 unique peptides (see Supplementary Table 1).Quantifying metagenomic presence of MAGsA reference database containing all 107 MAGs’ contigs was created (bbmap.sh command, default settings). Then, the paired-end short reads from each sample (FASTQ files) were mapped into the reference database (bbwrap.sh, default settings), producing alignment (SAM) files, which were converted into BAM format. Subsequently, a contig depth (coverage) table was produced using the command jgi_summarize_bam_contig_depths --outputDepth depth.txt --pairedContigs paired.txt *.bam. As each of the MAGs span on more than one contig, MAG depth in each sample was calculated as contig length weighted by the average depth. Finally, to account for unequal sequencing depth, each MAG depth was normalized to the number of short sequencing reads within the given sample.Correlating metagenomic and proteomic structuresIn order to compare metagenomic and proteomic structures, we first calculated the mean coding gene abundance and mean production levels of each of the 1629 detected core proteins over all 12 cows. Both mean gene abundance and mean production level were translated into ranks using the R rank function. The produced proteins were ranked in descending order and the coding genes in the gene abundance vector were reordered accordingly. The two reordered ranked vectors then plotted using the R pheatmap function, and colored using the same color scale.Selection of proteins for downstream analysisAs our goal was to analyze plasticity in microbial protein production in varying environments, e.g., as a function of host state, only MAGs that were identified in all of the 12 proteomic samples were kept for further analysis. Consequently, only proteins that were identified in at least half of the proteomic samples (e.g., in at least six samples) were selected. This last step aimed to reduce spurious correlation results. These filtering steps retained 79 MAGs coding for a total of 1,629 measurable proteins.Feed efficiency state prediction and ordinationIn order to calculate the accuracy in predicting host feed efficiency state based on the different data layers available (16S rRNA (Supplementary Table 2), metagenomics, metaproteomics), the principal component analysis (PCA) axes for all the samples based on the microbial protein production profiles were calculated. Then, twelve cycles of model building and prediction were made. Each time, the two first PCs of each of five cows along with their phenotype (efficiency state) were used to build a Support Vector Machine (SVM) [R caret package] prediction model and one sample was left out. The model was then used to perform subsequent prediction of the left-out animal phenotype (feed efficiency) by feeding the model with that animal’s first two PCs. This leave-one-out methodology was then repeated over all the samples. Finally, the prediction accuracy was determined as the percent of the cases where the correct label was assigned to the left-out sample. For the proteomics data, this procedure was applied on both the raw protein counts, and the protein production normalized based on MAG abundance, which enabled us to compare the prediction accuracies of the microbial protein production to that of the raw protein counts.Identification proteins associated with a specific host stateIn order to split the proteomics dataset into microbial proteins that tend to be produced differently as a function of the host feed efficiency states, each microbial protein profile was correlated to the sample’s host feed efficiency measure (as calculated by RFI) using the Spearman correlation (R function cor), disregarding the p value. Proteins that had a positive correlation to RFI were grouped as inefficiency associated proteins. In contrast, proteins that presented a negative correlation to RFI were grouped as efficiency associated proteins. To test for equal sizes of these two protein groups, a binomial test was performed (R function binom.test) to examine the probability to get a low number of feed efficient proteins from the overall proteins under examination, when the expected probability was set to 0.5.Functional assignment of proteinsProtein functions were assigned based on the KEGG (Kegg Encyclopedia of Genes and Genomes) [31] database. The entire KEGG genes database was compiled into a Diamond [32] search library. Then, the selected microbial proteins were searched against the database using the Diamond search tool. Significant hits (evalue < 5e-5) were further analyzed to identify the corresponding KO (KEGG Ortholog number). Annotations of glycoside hydrolases were performed using dbcan2 [33].Protein level checkerboard distribution across the feed efficiency groupsThe checkerboard distribution in protein production profiles was estimated separately within the feed efficient and inefficient animal groups. To enable the comparison between the two groups’ checkerboardness level, we chose a standardized C-score estimate (Standardized Effect Size C-score - S.E.S C-Score), based on the comparison of the observed C-score to a null-model distribution derived from simulations. The S.E.S C-score was estimated using the oecosimu function from R vegan package with 100,000 simulated null-model communities.Calculating functional redundancyThe functional redundancy within a given group of proteins was measured as the mean number of times a given KO occurred within a given group, while neglecting proteins that have not been assigned a KO level functional annotation.In order to test whether a given group of proteins exhibits more or less functional redundancy than would have been expected, a null distribution for functional redundancy was created, based on the number of proteins in the given group. A random group of proteins was drawn from the entire set, keeping the same sample size as in the tested group, and the process was repeated 100 times. Then, the functional redundancy for each random protein group was calculated. Thereafter, the null distribution was used to obtain a p value to measure the likelihood of obtaining such a value under the null.Examining functional divergenceExamining the functional divergence between the two groups of proteins, e.g. the feed efficiency and inefficiency associated proteins, was done by first counting the amount of shared functional annotations, in terms of KOs between the two groups. Thereafter, a null distribution for the expected count of KOs was built by randomly splitting in an iterative manner the proteins into groups of the same sizes and calculating the number of shared KOs. A p value for the actual count of shared proteins was obtained by ranking the actual count over the null distribution.Calculating average nearest neighbor ratio (ANN ratio)ANN Ratio analysis was carried out independently for each protein function (KO), containing more than 14 proteins with at least 5 proteins within each feed efficiency group. Initially, all proteins assigned to a given KO were split into two sets, in accordance to their feed efficiency affiliation group. Thereafter, proteins within each set were independently projected into two-dimensional space by PCA applied directly to Sequence Matrix [34]. Average nearest neighbor ratio within each set was then calculated within the minimum enclosing rectangle defined by principal component axes PC1 and PC2, as defined by Clark and Evans [35].MAG feed efficiency score calculationMicroorganism feed efficiency score was calculated for each MAG individually by first ranking each protein being produced by the given microbe along the 12 animals, based on the normalized protein production levels. Thereafter, a representative production value for the microbe in each animal was calculated as the average of the ranked (normalized) protein production levels in that animal (using R rank function). This ranking allowed us to alleviate the potential skewing effect of highly expressed proteins. The microorganism’s Feed Efficiency Score was calculated as the difference between its mean representative production value within feed efficient animals to that within feed inefficient animals. Values close to zero will reflect similar distribution between the two animal groups, positive values will indicate higher expression among efficient animals, and negative values will indicate higher expression among inefficient animals. To calculate significance, the actual feed efficiency score was compared to values in a distribution derived from a permutation based null model. Each of the permuted Feed Efficiency Scores (10,000 for each microbe) was obtained by independently shuffling each of the proteins produced by the MAG between the animals, prior to calculating the actual microorganism feed efficiency score. By positioning the absolute score value over its distribution under permuted assumptions (absolute values), we obtained a significance p value.MAG phylogenetic tree construction and phylogenetic signal estimationIn order to assess the link between phylogenetic similarity between the MAGs and their association with feed efficiency, phylogenetic tree estimating evolutionary relationships between the MAGs was constructed using the PhyloPhlAn pipeline [36]. The phylogenetic signal for Microorganism Feed Efficiency Score was estimated by providing the phylogSignal function from R phylosignal [37] package with MAGs phylogenetic tree and respective values. Pagel’s Lambda statistics was chosen for the analysis, owing to its robustness [38].Plot generationAll bar plots, scatter plots and other point plots were generated with R package ggplot2. Heatmaps were produced by either ggplot2 [39] or pheatmap [https://cran.r-project.org/web/packages/pheatmap/index.html] R packages. KEGG map was produced using the online KEGG Mapper tool [40]. Phylocorrelogram was produced with phyloCorrelogram function from R package phylosignal [37].MAG differential production analysisMAGs that contain a minimal number of proteins (50 functions) were selected for differential protein production analysis, in order to have sufficient data to perform statistical tests. For each MAG, the relative production was used in order to calculate the Jaccard pairwise dissimilarity for core protein production between feed efficient and inefficient cows using the R vegan package. Analysis of similarity between efficiency and inefficiency associated proteins for each MAG (ANOSIM) values and p values were then calculated using the same package.Predicting animal feed efficiency state according to GH family countsUsing all GH annotated proteins, a feature table that sums the count of each GH family within each sample was produced. Thereafter a leave-one-out cross-validation (LOOCV) [R caret package] was performed, each time building a Random Forest (RF) prediction model from the GH family counts and efficiency state of 11 samples, leaving one sample outside. Each one of the RF models, in its turn, was applied on the left-out animal to predict its efficiency state. Model accuracy and AUC curve were calculated based on the LOOCV performance. More

Azam, F. & Malfatti, F. Microbial structuring of marine ecosystems. Nat. Rev. Microbiol. 5, 782–791 (2007).CAS 
PubMed 
Article 

Google Scholar 
Martiny, J. B. H. et al. Microbial biogeography: Putting microorganisms on the map. Nat. Rev. Microbiol. 4, 102–112 (2006).CAS 
PubMed 
Article 

Google Scholar 
Hanson, C. A., Fuhrman, J. A., Horner-Devine, M. C. & Martiny, J. B. H. Beyond biogeographic patterns: Processes shaping the microbial landscape. Nat. Rev. Microbiol. 10, 497–506 (2012).CAS 
PubMed 
Article 

Google Scholar 
Grossart, H. P. Ecological consequences of bacterioplankton lifestyles: Changes in concepts are needed. Environ. Microbiol. Rep. 2, 706–714 (2010).PubMed 
Article 

Google Scholar 
Simon, M., Grossart, H. P., Schweitzer, B. & Ploug, H. Microbial ecology of organic aggregates in aquatic ecosystems. Aquat. Microb. Ecol. 28, 175–211 (2002).Article 

Google Scholar 
Smith, D. C., Simon, M., Alldredge, A. L. & Azam, F. Intense hydrolytic enzyme activity on marine aggregates and implication for rapid particle dissolution. Nature 359, 139–141 (1992).ADS 
CAS 
Article 

Google Scholar 
Grossart, H. P., Tang, K. W., Kiørboe, T. & Ploug, H. Comparison of cell-specific activity between free-living and attached bacteria using isolates and natural assemblages. FEMS Microbiol. Lett. 206, 194–200 (2007).Article 
CAS 

Google Scholar 
Rieck, A., Herlemann, D. P. R., Jürgens, K. & Grossart, H. Particle-associated differ from free-living bacteria in surface waters of the Baltic Sea. Front. Microbiol. 6, 1297 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
Karner, M. & Herndl, G. J. Extracellular enzymatic activity and secondary production in free-living and marine-snow-associated bacteria. Mar. Biol. 113, 341–347 (1992).CAS 
Article 

Google Scholar 
Lyons, M. M. & Dobbs, F. C. Differential utilization of carbon substrates by aggregate-associated and water-associated heterotrophic bacterial communities. Hydrobiologia 686, 181–193 (2012).CAS 
Article 

Google Scholar 
Simon, H. M., Smith, M. W. & Herfort, L. Metagenomic insights into particles and their associated microbiota in a coastal margin ecosystem. Front. Microbiol. 5, 466 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Smith, M. W., Allen, L. Z., Allen, A. E., Herfort, L. & Simon, H. M. Contrasting genomic properties of free-living and particle-attached microbial assemblages within a coastal ecosystem. Front. Microbiol. 4, 120 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
Mestre, M. et al. Spatial variability of marine bacterial and archaeal communities along the particulate matter continuum. Mol. Ecol. 26, 6827–6840 (2017).CAS 
PubMed 
Article 

Google Scholar 
Bižic-Ionescu, M. et al. Comparison of bacterial communities on limnic versus coastal marine particles reveals profound differences in colonization. Environ. Microbiol. 17, 3500–3514 (2015).PubMed 
Article 
CAS 

Google Scholar 
Hollibaugh, J. T., Wong, P. S. & Murrell, M. C. Similarity of particle-associated and free-living bacterial communities in northern San Francisco Bay, California. Aquat. Microb. Ecol. 21, 103–114 (2000).Article 

Google Scholar 
Ortega-Retuerta, E., Joux, F., Jeffrey, W. H. & Ghiglione, J. F. Spatial variability of particle-attached and free-living bacterial diversity in surface waters from the Mackenzie River to the Beaufort Sea (Canadian Arctic). Biogeosciences 10, 2747–2759 (2013).ADS 
Article 

Google Scholar 
Noble, P. A., Bidle, K. D. & Fletcher, M. Natural microbial community compositions compared by a back-propagating neural network and cluster analysis of 5S rRNA. Appl. Environ. Microbiol. 63, 1762–1770 (1997).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Zhou, J. & Ning, D. Stochastic community assembly: Does it matter in microbial ecology?. Microbiol. Mol. Biol. Rev. 81, e00002-17 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Jain, A., Balmonte, J. P., Singh, R., Bhaskar, P. V. & Krishnan, K. P. Spatially resolved assembly, connectivity and structure of particle-associated and free-living bacterial communities in a high Arctic fjord. FEMS Microbiol. Ecol. 97, 1–12 (2021).Article 
CAS 

Google Scholar 
Yao, Z. et al. Bacterial community assembly in a typical estuarine marsh. Appl. Environ. Microbiol. 85, e02602-18 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Wang, J. et al. Assembly processes and source tracking of planktonic and benthic bacterial communities in the Yellow River estuary. Environ. Microbiol. 23, 2578–2591 (2021).PubMed 
Article 

Google Scholar 
Balmonte, J. P. et al. Sharp contrasts between freshwater and marine microbial enzymatic capabilities, community composition, and DOM pools in a NE Greenland fjord. Limnol. Oceanogr. 65, 77–95 (2019).ADS 
Article 
CAS 

Google Scholar 
Fortunato, C. S., Herfort, L., Zuber, P., Baptista, A. M. & Crump, B. C. Spatial variability overwhelms seasonal patterns in bacterioplankton communities across a river to ocean gradient. ISME J. 6, 554–563 (2012).CAS 
PubMed 
Article 

Google Scholar 
Yawata, Y., Carrara, F., Menolascina, F. & Stocker, R. Constrained optimal foraging by marine bacterioplankton on particulate organic matter. Proc. Natl. Acad. Sci. USA 117, 25571–25579 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hu, Y. et al. The relationships between the free-living and particle-attached bacterial communities in response to elevated eutrophication. Front. Microbiol. 11, 423 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Lima-Mendez, G. et al. Determinants of community structure in the grobal plankton interactome. Science (80-) 348, 1262073-1–10 (2015).Article 
CAS 

Google Scholar 
Milici, M. et al. Co-occurrence analysis of microbial taxa in the Atlantic ocean reveals high connectivity in the free-living bacterioplankton. Front. Microbiol. 7, 649 (2016).PubMed 
PubMed Central 

Google Scholar 
Herren, C. M. & McMahon, K. D. Cohesion: A method for quantifying the connectivity of microbial communities. ISME J. 11, 2426–2438 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Friedman, J. & Alm, E. J. Inferring correlation networks from genomic survey data. PLoS Comput. Biol. 8, e1002687 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Deng, Y. et al. Molecular ecological network analyses. BMC Bioinform. 13, 113 (2012).Article 

Google Scholar 
Labry, C. et al. High alkaline phosphatase activity in phosphate replete waters: The case of two macrotidal estuaries. Limnol. Oceanogr. 61, 1513–1529 (2016).ADS 
Article 

Google Scholar 
Crump, B. C. et al. Quantity and quality of particulate organic matter controls bacterial production in the Columbia River estuary. Limnol. Oceanogr. 62, 2713–2731 (2017).ADS 
CAS 
Article 

Google Scholar 
Canuel, E. A. & Hardison, A. K. Sources, ages, and alteration of organic matter in Estuaries. Ann. Rev. Mar. Sci. 8, 409–434 (2016).PubMed 
Article 

Google Scholar 
He, W., Chen, M., Schlautman, M. A. & Hur, J. Dynamic exchanges between DOM and POM pools in coastal and inland aquatic ecosystems: A review. Sci. Total Environ. 551–552, 415–428 (2016).ADS 
PubMed 
Article 
CAS 

Google Scholar 
Bianchi, T. S. The role of terrestrially derived organic carbon in the coastal ocean: A changing paradigm and the priming effect. Proc. Natl. Acad. Sci. 108, 19473–19481 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Auffret, G. A. Dynamique sédimentaire de la marge continentale celtique-Evolution Cénozoïque-Spécificité du Pleistocène supérieur et de l’Holocène (Université de Bordeaux I, 1983).
Google Scholar 
Delmas, R. & Tréguer, P. Évolution saisonnière des nutriments dans un écosystème eutrophe d’Europe occidentale (la rade de Brest). Interactions marines et terrestres. Oceanol. Acta 6, 345–356 (1983).CAS 

Google Scholar 
Bassoullet, P. Etude de la dynamique des sédiments en suspension dans l’estuaire de l’Aulne (rade de Brest) (Université de Bretagne Occidentale, 1979).
Google Scholar 
Klindworth, A. et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 41, e1 (2013).CAS 
PubMed 
Article 

Google Scholar 
Bolyen, E. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Olesen, S. W., Duvallet, C. & Alm, E. J. dbOTU3: A new implementation of distribution-based OTU calling. PLoS ONE 12, 1–13 (2017).Article 
CAS 

Google Scholar 
Quast, C. et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41, 590–596 (2013).Article 
CAS 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (2013).Whickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).Book 

Google Scholar 
Lê, S., Josse, J. & Husson, F. FactoMineR: An R package for multivariate analysis. J. Stat. Softw. 25, 1–18 (2008).Article 

Google Scholar 
Wei, T. & Simko, V. R package ‘corrplot’: Visualization of a Correlation Matrix (2011).McMurdie, P. J. & Holmes, S. Phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8, e61217 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Oksanen, J. et al. Vegan: Community Ecology Package (2022).Liu, C., Cui, Y., Li, X. & Yao, M. Microeco: An R package for data mining in microbial community ecology. FEMS Microbiol. Ecol. 97, fiaa255 (2021).CAS 
PubMed 
Article 

Google Scholar 
Kandlikar, G. ranacapa: Utility Functions and ‘shiny’ App for Simple Environmental DNA Visualizations and Analyses (2021).Cao, Y. microbiomeMarker: microbiome biomarker analysis toolkit (2021).Tsirogiannis, C. & Brody, S. PhyloMeasures: Fast and Exact Algorithms for Computing Phylogenetic Biodiversity Measures (2017).McKnight, D. T. et al. Methods for normalizing microbiome data: An ecological perspective. Methods Ecol. Evol. 10, 389–400 (2019).Article 

Google Scholar 
Paradis, E. & Schliep, K. Ape 50: An environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics https://doi.org/10.1093/bioinformatics/bty633 (2019).Article 
PubMed 

Google Scholar 
Legendre, P. & Legendre, L. Numerical Ecology (Third English Edition) (Elsevier, 2012).MATH 

Google Scholar 
Stegen, J. C., Lin, X., Fredrickson, J. K. & Konopka, A. E. Estimating and mapping ecological processes influencing microbial community assembly. Front. Microbiol. 6, 370 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
Stegen, J. C., Lin, X., Konopka, A. E. & Fredrickson, J. K. Stochastic and deterministic assembly processes in subsurface microbial communities. ISME J. 6, 1653–1664 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Naimi, B. usdm: Uncertainty Analysis for Species Distribution Models (2017).Wu, W., Xu, Z., Dai, M., Gan, J. & Liu, H. Homogeneous selection shapes free-living and particle-associated bacterial communities in subtropical coastal waters. Divers. Distrib. 00, 1–14 (2020).
Google Scholar 
Wang, Y. et al. Patterns and processes of free-living and particle-associated bacterioplankton and archaeaplankton communities in a subtropical river-bay system in South China. Limnol. Oceanogr. 65, 161–179 (2020).
Google Scholar 
Zhou, L. et al. Environmental filtering dominates bacterioplankton community assembly in a highly urbanized estuarine ecosystem. Environ. Res. 196, 110934 (2021).CAS 
PubMed 
Article 

Google Scholar 
Graham, E. B. & Stegen, J. C. Dispersal-based microbial community assembly decreases biogeochemical function. Processes 5, 65 (2017).Article 

Google Scholar 
Campbell, B. J. & Kirchman, D. L. Bacterial diversity, community structure and potential growth rates along an estuarine salinity gradient. ISME J. 7, 210–220 (2013).CAS 
PubMed 
Article 

Google Scholar 
Herlemann, D. P. R. et al. Transitions in bacterial communities along the 2000 km salinity gradient of the Baltic Sea. ISME J. 5, 1571–1579 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Fuhrman, J. A., Cram, J. A. & Needham, D. M. Marine microbial community dynamics and their ecological interpretation. Nat. Rev. Microbiol. 13, 133–146 (2015).CAS 
PubMed 
Article 

Google Scholar 
Sunagawa, S. et al. Structure and function of the global ocean microbiome. Science (80-) 348, 1261359 (2015).Article 
CAS 

Google Scholar 
Buchan, A., LeCleir, G. R., Gulvik, C. A. & González, J. M. Master recyclers: features and functions of bacteria associated with phytoplankton blooms. Nat. Rev. Microbiol. 12, 686–698 (2014).CAS 
PubMed 
Article 

Google Scholar 
Martinez-Garcia, M. et al. Capturing single cell genomes of active polysaccharide degraders: An unexpected contribution of verrucomicrobia. PLoS ONE 7, e35314 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Reintjes, G., Arnosti, C., Fuchs, B. M. & Amann, R. An alternative polysaccharide uptake mechanism of marine bacteria. ISME J. 11, 1640–1650 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Gralka, M., Szabo, R., Stocker, R. & Cordero, O. X. Trophic interactions and the drivers of microbial community assembly. Curr. Biol. 30, R1176–R1188 (2020).CAS 
PubMed 
Article 

Google Scholar 
Liu, J., Meng, Z., Liu, X. & Zhang, X. H. Microbial assembly, interaction, functioning, activity and diversification: a review derived from community compositional data. Mar. Life Sci. Technol. 1, 112–128 (2019).ADS 
Article 

Google Scholar 
Hernandez, D. J., David, A. S., Menges, E. S., Searcy, C. A. & Afkhami, M. E. Environmental stress destabilizes microbial networks. ISME J. 15, 1722–1734 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Herren, C. M. & McMahon, K. D. Keystone taxa predict compositional change in microbial communities. Environ. Microbiol. 20, 2207–2217 (2018).PubMed 
Article 

Google Scholar 
Liénart, C. et al. Dynamics of particulate organic matter composition in coastal systems: A spatio-temporal study at multi-systems scale. Prog. Oceanogr. 156, 221–239 (2017).Article 

Google Scholar 
Fraisse, S., Bormans, M. & Lagadeuc, Y. Morphofunctional traits reflect differences in phytoplankton community between rivers of contrasting flow regime. Aquat. Ecol. 47, 315–327 (2013).Article 

Google Scholar 
Treguer, P. & Queguiner, B. Seasonal variations in conservative and nonconservative mixing of nitrogen compounds in a West European macrotidal estuary. Oceanol. Acta 12, 371–380 (1989).CAS 

Google Scholar 
Grossart, H. P. & Tang, K. W. Communicative & integrative biology. Commun. Integr. Biol. 3, 491–494 (2010).PubMed 
PubMed Central 
Article 

Google Scholar  More

Taxonomy, variation, and biochronologyThe fossils described herein represent one of the most valuable and best-preserved samples of “Dama-like” deer from the European Early Pleistocene. The systematics of these forms has been essentially based on the morphology of the antlers and teeth, with less attention paid to the skull (due to the rarity of well-preserved finds) and postcranial bones.The Pantalla sample shows a combination of characters allowing an unambiguous attribution to ‘P.’ nestii, a species reported confidently so far in the early Late Villafranchian of Italy (several sites) and in the Georgian Homo-bearing locality of Dmanisi (Supplementary Table S1). Based on the literature6,8,12,38, these characters include: four-pointed antler with elongated, slender, and tubular beam; basal tine branching off at a certain distance from the burr forming an acute angle; well-developed middle tine; terminal bifurcation oriented normal to the sagittal plane; cranium with large orbits, preorbital fossae, and ethmoidal vacuities; relatively elongated neurocranium with flat parietals; caudally-oriented pedicles; molarized P2-P3; presence of cingula in upper molars; enlarged i1; un-molarized p4. However, some characters observed in the Pantalla specimens (e.g., rostral edge of the orbit reaching the level of M2; elongated metapodials) do not fit the revised diagnosis of ‘P.’ nestii by Croitor12. The latter author considers nestii as the earliest species of the genus Cervus based on similarities with the extant red deer especially in cranial morphology12,22,23. However, in our opinion, his conclusions are biased by relying mostly on the skull IGF 243 of ‘P.’ nestii from Upper Valdarno6,8, which is heavily deformed and belongs to a juvenile individual (see below for details on ontogenetic variation in ‘Pseudodama’).A broader look at the entire record of ‘P.’ nestii reveals that this species displays a mosaic of characters between Dama and Cervus, but also that the shared characters with Dama are prevalent (as already pointed out by Azzaroli8). The Pantalla sample allows to substantiate these conclusions very well. Our CT-based comparisons between the crania from Pantalla and those of extant red deer and fallow deer (Fig. 3) highlight some morphological similarities with the former, including a relatively longish neurocranium with steep forehead and deep preorbital fossa. On the other hand, ‘P.’ nestii from Pantalla clearly shows Dama-like cranial characters, such as a marked interfrontal crest, horizontal zygomatic arch, high maxilla below the orbit, muzzle more inclined ventrally and less cylindrical in overall shape, sub-horizontal upper cheek tooth row (i.e., the occlusal margin of the row is approximately straight in buccal view), apical surface of the pedicle more inclined dorsocaudally, and overall morphology of the antlers, which in rostral view diverge, rather than converge as in the red deer (Fig. 2).Likewise, the teeth from Pantalla, have a mixture of Dama and Cervus characters although the former are prevalent. All the premolar characters (the complete absence of a lingual grove on P4, the presence of a cingulum on the distolingual wall of P4, the presence of a small paraconid in p2, the entoconid more aligned with the mesiodistal axis in p3-p4, and a weak mesial cingulum on p4) and most of the lower molar characters are Dama-like. The upper molar features are instead more reminiscent of Cervus being either intermediate between the morphology of the latter and that of extant Dama or even matching Cervus (see Supplementary Table S6 and below).The postcranial remains from Pantalla appear more similar to Dama than to Cervus. Of the 23 morphological characters by Lister39 which are present in the preserved bones (axis, metacarpal, tibia, astragalus, calcaneum, cubo-navicular, metatarsal, phalanx I, and phalanx II), 21 scores as fallow deer and only two as red deer (details in Supplementary Table S7).A mixed character suite between Dama and Cervus are revealed also by our palaeoneurological analysis. The brain of ‘P.’ nestii shows Dama-like size and Cervus-like morphology with a prominent cerebellum and a dorsoventrally flattened cerebrum. The latter character is clearly noticeable in ‘Pseudodama’ and Eucladoceros, is less evident in extant Cervus, and is missing in Dama. The hypothesis that depressed and longish cerebra represent a primitive character in Cervini (at least in Pleistocene European forms) is supported by our preliminary data and agree with Azzaroli8.Most interestingly, the two crania from Pantalla actually show some remarkable morphological differences. The neurocranium of 337643 is more lengthened (i.e., more Cervus-like), albeit this shape might be taphonomically modified by the lateral compression of the specimen. This morphology fits that observed in some other ‘P.’ nestii specimens such as IGF 1403 from Olivola (Italy), while the relatively shorter and more rounded neurocranium of 337655 resembles that of other specimens such as IGF 1404 also from Olivola. Moreover, 337643 shows a stronger nuchal crest than 337655. These differences may be related to ontogenesis (see the advanced age of 337643 based on tooth wear). In several cervid species including fallow deer, aging leads to morphological changes in the neurocranium, which tends to elongate and flatten and shows a more developed nuchal crest, probably as a response to the support of larger and heavier antlers18,38. Similarly, in 337643, the pedicles are apparently closer to one another due to their thickening—an expected condition for an old individual as the distance between the pedicles tends to decrease with age8—and markedly shorter than wide. Our comparative data on European Dama-like deer show that the pedicle section can be highly variable both within and between species, although a general trend of laterolateral flattening (i.e., oval shape with major axis oriented anteroposteriorly) can be traced through time (Supplementary Fig. S3), probably as a result of the development of wide, laterally-projecting palmated antlers (in extant deer, D. dama is among those with the heaviest antlers relative to body size40,41). Therefore, the Pantalla sample on the one hand confirms the variation in cranial morphology already observed for ‘P.’ nestii6,8, on the other hand it supports the affinities between this species and the fallow deer. The presence of Cervus-like features especially in cranial morphology may be interpreted as plesiomorphic characters which, associated with some characters of the dentition and of the brain, suggest a basal position of ‘Pseudodama’ in the evolutionary history of the Cervini. This hypothesis may be tested in the future through phylogenetic analyses, currently made difficult by the lack of sufficiently well-preserved material of some species of ‘Pseudodama’ (e.g., ‘P.’ lyra, ‘P.’ perolensis).Compared with other specimens of ‘P.’ nestii6,8, the sample from Pantalla shows some plesiomorphic characters including a high ratio between the premolar and molar lengths, i.e., 0.77–0.82 (n = 2) for upper teeth (LP/LM) and 0.68–0.69 (n = 3) for lower teeth (Lp/Lm). These values are closer to the basal forms of ‘Pseudodama’, such as ‘P.’ lyra from Montopoli (LP/LM = 0.73, n = 1; Lp/Lm = 0.64, n = 2) and ‘P.’ rhenana from Saint Vallier (LP/LM = 0.75, n = 9; Lp/Lm = 0.68, n = 18; data from Valli42), than to ‘P.’ nestii from Olivola and Upper Valdarno (LP/LM = 0.72, n = 10; Lp/Lm = 0.63, n = 17). Other putatively plesiomorphic features of the sample from Pantalla are all those that approach it morphologically to Cervus (see Supplementary Table S6), i.e., the strong development of lingual conids and stylids in lower molars (Char. 439) and of buccal cones and styles in upper molars (Char. 139), the lack of a clear step between 2nd and 3rd lobe of m3 (Char. 1139), the strong lingual cingulum on upper molars (Char. 339), and the lack of the horizontal turning of the buccal columns of upper molars (Char. 439—the so-called buccal “cingulum”43). The strong lingual cingulum on upper molars is constantly present in the earliest species of the ‘Pseudodama’ group, ‘P.’ pardinensis9, and still present, although extremely rare, in ‘P.’ lyra from Montopoli, ‘P.’ rhenana from Saint Vallier and Senèze, ‘P.’ perolensis from Peyrolles, and ‘P.’ nestii from Olivola. However, this feature is back less rare in ‘P.’ nestii from Upper Valdarno and ‘P.’ farnetensis from Selvella, suggesting a certain polymorphism at this stage. The lack of buccal “cingulum” is a constant in the earliest ‘Pseudodama’ populations (Lower Valdarno, Saint Vallier, Senèze), the buccal “cingulum” appearing, although rare, in ‘P.’ perolensis from Peyrolles and ‘P.’ nestii from Olivola and Upper Valdarno but becoming more common only in later ‘P.’ farnetensis, ‘P.’ vallonnetensis, and constant in Dama.The above affinities between the Pantalla deer and the early representatives of ‘Pseudodama’ support the idea that the age of the assemblage may be close to the beginning of the Late Villafranchian (ca. 2.1–2.0 Ma), as already suggested based on the occurrence of Leptobos merlai44 and a primitive form of Equus stenonis35. Thus, the ‘P.’ nestii sample described herein may represent one of the earliest occurrences of the species in Europe.Palaeoecological and palaeoethological inferencesThe Pantalla sample is also noteworthy as it allows opening a window into the behaviour of these extinct deer. The anomalies found on the two male crania are probably the result of different traumas during their life.Deer are well known for the intense fights they engage in during the rutting season using their antlers, as a result of an escalation of a broad repertoire of threats and displays45. Mineralized antlers are solid structures able to withstand the vehemence of the fight46, whereas growing antlers are extremely fragile and any contact with a solid object may result in a serious injury47,48 that may jeopardize the bearer’s ability to compete with conspecifics and, consequently, its dominance status49. Accidents are inevitable in the life of a deer and, in case of the suffered damage not leading to the breakage of the growing beam and consequent loss of its distal part, the antler may continue its growth although, in case of a severe lesion, at a crooked angle45. Thus, if the antler was just cracked and the broken part was held together by the velvet and periosteum, with the blood supply still being guaranteed, the damaged beam would just present a conspicuous swelling around the area of fracture (i.e., a fracture callus)45,50 and a change in the axis of orientation. These features match those seen in the left beam of 337655, which shows a fracture callus between the basal and middle tines corresponding to a change in the orientation of the beam.The supernumerary tine of the right antler of the same individual can be interpreted as the result of a trauma, too. Considering the delicate nature of the growing antlers and the non-negligible risks of occurrence of an injury, it is safe to believe that the right antler has undergone a light traumatic event (most likely concerning the pedicle) at some early stage of its growth. In fact, it is known that limited injuries could result in the growth of supernumerary tines, even in atypical positions51, as it has been documented in other deer species (e.g., reindeer52, sambar53). It is therefore reasonable to hypothesize that both antler anomalies of 337655 derive from traumas suffered by the deer during the antler growth, when the velvet was still present. It is not known whether the two injuries happened at the same time or in two different events. In fact, it cannot even be said that the two events took place during the same season. While the breakage of the left beam must have occurred in the year of the animal’s death (i.e., during the velvet period preceding the period of hard antler in which the individual died), the development of the supernumerary tine on the right may be the result of a trauma suffered in a previous year. This is due to the fact that when unilateral trauma affects the generative region of the antler (i.e., the pedicle area), abnormalities such as supernumerary tines can reappear in next antler cycles even in more intensified forms54, as in the case of 337655 in which the extra-tine is extremely long.The bone anomaly on the right squamosal of 337643 is also likely the outcome of an injury. Although the external portion was artificially smoothed during the preparation of the specimen, the outer and inner morphology matches that of a callus related to the healing of a major lesion and probable intracranial abscessation. Post-traumatic inflammatory processes are known to cause erosion or pitting of cranial bones in deer55 and can be triggered by many factors (e.g., wounds and abrasions of the pedicle56), among which violent sexual competition among males with hard antlers is considered one of the most common55,57. The advanced healing of the injury shown by 337643 suggests that it was not the cause of death, but rather that the individual survived a long time after the trauma albeit with the brain partially compressed by the callus.The six mandibles recovered at Pantalla, all coming from the same bone accumulation hence reasonably referable to a single deer population, represent several age classes, from calves as young as a few months up to very old individuals (i.e., over 15 years; Supplementary Table S4). Unfortunately, no mandible can be safely associated with the two male crania, although 337631 may belong to the same individual as 337643 based on advanced wear and size. Interestingly, the three most significant cranial remains (crania 337655 and 337643 and frontal bone fragment with basal antler base 337625) belong to adult males, which probably died during the hard antler period (i.e., rutting season: 337655 and 337625) or shortly after (i.e., 337643). The absence of females (at least among the remains with certain sex attribution) contrasts with the population structure in the extant fallow deer, in which females represent on average 75% of the herd58. However, the relative abundance of males may increase up to 50% in the rutting season59,60. Therefore, in spite of the relatively low number of fossils available, based on the age and sex structure of the palaeopopulation and by analogy with the extant fallow deer, the most plausible hypothesis is that the Pantalla deer died during or immediately after the rutting season (Fig. 5).Figure 5Life appearance of ‘Pseudodama’ nestii represented during the rutting season. The reconstruction is based on the cranial and postcranial material from the Early Pleistocene of Pantalla (Italy) and on literature data. Artwork by D.A. Iurino.Full size image More

Wiens, J. J. Climate-related local extinctions are already widespread among plant and animal species. PLoS Biol. 14, e2001104 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Scheffers, B. R. et al. The broad footprint of climate change from genes to biomes to people. Science 354, aaf7671 (2016).PubMed 
Article 
CAS 

Google Scholar 
Parmesan, C. & Yohe, G. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42 (2003).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Urban, M. C. Accelerating extinction risk from climate change. Science 348, 571–573 (2015).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Wingfield, J. C. et al. Organism-environment interactions in a changing world: a mechanistic approach. J. Ornithol. 152, 279–288 (2011).Article 

Google Scholar 
Mendoza-Gonzalez, G., Martinez, M. L., Rojas-Soto, O. R., Vazquez, G. & Gallego-Fernandez, J. B. Ecological niche modeling of coastal dune plants and future potential distribution in response to climate change and sea level rise. Glob. Change Biol. 19, 2524–2535 (2013).ADS 
Article 

Google Scholar 
Saunders, S. P. et al. Community science validates climate suitability projections from ecological niche modeling. Ecol. Appl. 30, 17 (2020).Article 

Google Scholar 
Peterson, A. T., Cobos, M. E. & Jimenez-Garcia, D. Major challenges for correlational ecological niche model projections to future climate conditions. Ann. N. Y. Acad. Sci. 1429, 66–77 (2018).ADS 
PubMed 
Article 

Google Scholar 
Mays, H. L. et al. Genomic analysis of demographic history and Ecological niche modeling in the endangered Sumatran Rhinoceros Dicerorhinus sumatrensis. Curr. Biol. 28, 70–76 (2018).CAS 
PubMed 
Article 

Google Scholar 
Malcolm, R. J., Liu, C., Neilson, P. R., Hansen, L. & Hannah, L. A. Global warming and extinctions of endemic species from biodiversity hotspots. Conserv. Biol. 20, 538–548 (2005).Article 

Google Scholar 
Fitzpatrick, M. C. & Keller, S. R. Ecological genomics meets community-level modelling of biodiversity: mapping the genomic landscape of current and future environmental adaptation. Ecol. Lett. 18, 1–16 (2015).PubMed 
Article 

Google Scholar 
Gotelli, J. N. & Stanton-Geddes, J. Climate change, genetic markers and species distribution modelling. J. Biogeogr. 42, 1577–1585 (2015).Article 

Google Scholar 
Ruegg, K. et al. Ecological genomics predicts climate vulnerability in an endangered southwestern songbird. Ecol. Lett. 21, 1085–1096 (2018).PubMed 
Article 

Google Scholar 
Razgour, O. et al. Considering adaptive genetic variation in climate change vulnerability assessment reduces species range loss projections. Proc. Natl Acad. Sci. USA 116, 10418–10423 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Valladares, F. et al. The effects of phenotypic plasticity and local adaptation on forecasts of species range shifts under climate change. Ecol. Lett. 17, 1351–1364 (2014).PubMed 
Article 

Google Scholar 
Bay, R. A. et al. Genomic signals of selection predict climate-driven population declines in a migratory bird. Science 359, 83–86 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Rhone, B. et al. Pearl millet genomic vulnerability to climate change in West Africa highlights the need for regional collaboration. Nat. Commun. 11, 5274 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Rahbek, C. et al. Building mountain biodiversity: geological and evolutionary processes. Science 365, 1114–1119 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Fjeldså, J., Bowie, R. C. K. & Rahbek, C. The role of mountain ranges in the diversification of birds. Annu. Rev. Ecol. Evol. Syst. 43, 249–265 (2012).Article 

Google Scholar 
Freeman, B. G., Scholer, M. N., Ruiz-Gutierrez, V. & Fitzpatrick, J. W. Climate change causes upslope shifts and mountaintop extirpations in a tropical bird community. Proc. Natl Acad. Sci. USA 115, 11982–11987 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
He, J. K., Lin, S. L., Li, J. T., Yu, J. H. & Jiang, H. S. Evolutionary history of zoogeographical regions surrounding the Tibetan Plateau. Commun. Biol. 3, 9 (2020).Article 
CAS 

Google Scholar 
Wu, Y. J. et al. Explaining the species richness of birds along a subtropical elevational gradient in the Hengduan Mountains. J. Biogeogr. 40, 2310–2323 (2013).Article 

Google Scholar 
del Hoyo, J., Elliott, A., Sargatal, J. & Christie, D. A. Handbook of the Birds of the World (Lynx Edicions, 2013).Qu, Y. et al. Lineage diversification and historical demography of a montane bird Garrulax elliotii – implications for the Pleistocene evolutionary history of the eastern Himalayas. BMC Evolut. Biol. 11, 174 (2011).Article 

Google Scholar 
Qu, Y. et al. Long-term isolation and stability explain high genetic diversity in the Eastern Himalaya. Mol. Ecol. 23, 705–720 (2014).PubMed 
Article 

Google Scholar 
Wang, W. J. et al. Glacial expansion and diversification of an East Asian montane bird, the green-backed tit (Parus monticolus). J. Biogeogr. 40, 1156–1169 (2013).Article 

Google Scholar 
Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210–3212 (2015).PubMed 
Article 
CAS 

Google Scholar 
Laine, V. N. et al. Evolutionary signals of selection on cognition from the great tit genome and methylome. Nat. Commun. 7, 9 (2016).Article 
CAS 

Google Scholar 
Ellis, N., Smith, S. J. & Pitcher, C. R. Gradient forests: calculating importance gradients on physical predictors. Ecology 93, 156–168 (2012).PubMed 
Article 

Google Scholar 
Giorgetta, M. A. Climate and carbon cycle changes from 1850 to 2100 in MPI-ESM simulations for the Coupled Model Intercomparison Project phase 5. J. Adv. Model. Earth Syst. 5, 572–597 (2013).ADS 
Article 

Google Scholar 
Gent, P. R. et al. The community climate system model version 4. J. Clim. 24, 4973–4991 (2011).ADS 
Article 

Google Scholar 
Watanabe, M. et al. Improved climate simulation by MIROC5: mean states, variability, and climate sensitivity. J. Clim. 23, 6312–6335 (2010).ADS 
Article 

Google Scholar 
Voldoire, A. et al. The CNRM-CM5.1 global climate model: description and basic evaluation. Clim. Dyn. 40, 2091–2121 (2013).Article 

Google Scholar 
Frichot, E., Schoville, S. D., Bouchard, G. & Francois, O. Testing for associations between loci and environmental gradients using latent factor mixed models. Mol. Biol. Evol. 30, 1687–1699 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Forester, B. R., Jones, M. R., Joost, S., Landguth, E. L. & Lasky, J. R. Detecting spatial genetic signatures of local adaptation in heterogeneous landscapes. Mol. Ecololgy 25, 104–120 (2016).CAS 
Article 

Google Scholar 
Forester, B. R., Lasky, J. R., Wagner, H. H. & Urban, D. L. Comparing methods for detecting multilocus adaptation with multivariate genotype-environment associations. Mol. Ecol. 27, 2215–2233 (2018).CAS 
PubMed 
Article 

Google Scholar 
Li, C. et al. Two Antarctic penguin genomes reveal insights into their evolutionary history and molecular changes related to the Antarctic environment. Gigascience 3, 27 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Pirri, F. et al. Selection-driven adaptation to the extreme Antarctic environment in Emperor penguin. Preprint at bioRxiv https://doi.org/10.1101/2021.12.14.471946 (2021).Wang, L. C. et al. Involvement of the Arabidopsis HIT1/AtVPS53 tethering protein homologuein the acclimation of the plasma membrane to heat stess.J. Exp. Bot. 62, 3609–3620 (2011).CAS 
PubMed 
Article 

Google Scholar 
Piñol, R. A. et al. Preoptic BRS3 neurons increase body temperature and heart rate via multiple pathways. Cell Metab. 33, 1389–1403 (2021).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Guilherme, A. et al. Neuronal modulation of brown adipose activity through perturbation of white adipocyte lipogenesis. Mol. Metab. 16, 116–125 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Zhang, Y., Guo, W., zhang, Y., Zhang, H. & Wu, C. Insights into hypoxic adaptation in Tibetan chicken embryos from comparative proteomics. Comp. Biochem. Physiol. Part D. 31, 100602 (2019).CAS 

Google Scholar 
Pizzagalli, M. D., Bensimon, A. & Superti-Furga, G. A guide to plasma membrane solute carrier proteins. FEBS J. 288, 2784–2835 (2020).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Qu, Y. et al. Rapid phenotypic evolution with shallow genomic differentiation during early stages of high elevation adaptation in Eurasian Tree Sparrows. Natl Sci. Rev. 7, 113–127 (2020).PubMed 
Article 

Google Scholar 
Alexander, D. H., Novembre, J. & Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ferrier, S., Manion, G., Elith, J. & Richardson, K. Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment. Diversity Distrib. 13, 252–264 (2007).Article 

Google Scholar 
Thuiller, W., Lafourcade, B., Engler, R. & Araújo, M. B. BIOMOD – a platform for ensemble forecasting of species distributions. Ecography 32, 369–373 (2009).Article 

Google Scholar 
Chen, Y. et al. Large-scale genome-wide reveals climate adaptive variability in a cosmopolitan pest. Nat. Commun. 12, 7206 (2021).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Clarke, R. T., Rothery, P. & Raybould, A. F. Confidence limits for regression relationships between distance matrices: Estimating gene flow with distance. J. Agric. Biol. Environ. Stat. 7, 361–372 (2002).Article 

Google Scholar 
Excoffier, L., Dupanloup, I., Huerta-Sanchez, E., Sousa, V. C. & Foll, M. Robust demographic inference from genomic and SNP data. PLoS Genet. 9, e1003905 (2013).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Foden, W. B. et al. Climate change vulnerability assessment of species. WIREs Clim. Change 10, e551 (2019).Article 

Google Scholar 
Smith, T. B. et al. Genomic vulnerablity and soci-economic threats under climate change in an African rainforest bird. Evolut. Appl. 14, 1239–1247 (2021).Article 

Google Scholar 
Liu, B., Liang, E. Y., Liu, K. & Camarero, J. J. Species- and elevation-dependent growth responses to climate warming of mountain forests in the Qinling Mountains, central China. Forests 9, 11 (2018).
Google Scholar 
Dang, H. S., Zhang, Y. J., Zhang, K. R., Jiang, M. X. & Zhang, Q. F. Climate-growth relationships of subalpine fir (Abies fargesii) across the altitudinal range in the Shennongjia Mountains, central China. Clim. Change 117, 903–917 (2013).ADS 
Article 

Google Scholar 
Lingua, E., Cherubini, P., Motta, R. & Nola, P. Spatial structure along an altitudinal gradient in the Italian central Alps suggests competition and facilitation among coniferous species. J. Veg. Sci. 19, 425–436 (2008).Article 

Google Scholar 
Zhang, D. C., Zhang, Y. H., Boufford, D. E. & Sun, H. Elevational patterns of species richness and endemism for some important taxa in the Hengduan Mountains, southwestern China. Biodivers. Conserv. 18, 699–716 (2009).Article 

Google Scholar 
Zhang, R. Z., Zheng, D., Yang, Q. Y. & Liu, Y. H. Physical Geography of Hengduan Mountains (Science Press, 1997).Liu, Y. et al. Sino-Himalayan mountains act as cradles of diversity and immigration centres in the diversification of parrotbills (Paradoxornithidae). J. Biogeogr. 43, 1488–1501 (2016).Bush, A. et al. Incorporating evolutionary adaptation in species distribution modeling reduces projected vulnerability to climate change. Ecol. Lett. 17, 1468–148 (2016).Article 

Google Scholar 
Sparks, M. M., Westley, A. A. H., Falke, J. A. & Quinn, T. P. Thermal adaptation and phenotypic plasticity in a warming world: insights from common garden experiments on Alaskan sockeye salmon. Glob. Change Biol. 23, 5203–5217 (2017).ADS 
Article 

Google Scholar 
Merow, C., Wilson, A. M. & Jetz, W. Integrating occurrence data and expert maps for improved species range predictions. Glob. Ecol. Biogeogr. 26, 243–258 (2017).Article 

Google Scholar 
Weisenfeld, N. I., Kumar, V., Shah, P., Church, D. M. & Jaffe, D. B. Direct determination of diploid genome sequences. Genome Res. 27, 757–767 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Luo, R. et al. SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. Gigascience 1, 18 (2012).PubMed 
PubMed Central 
Article 

Google Scholar 
Altschul, S. F. et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 25, 3389–3402 (1997).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
She, R., Chu, J. S. C., Wang, K., Pei, J. & Chen, N. GenBlastA: enabling BLAST to identify homologous gene sequences. Genome Res. 19, 143–149 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Birney, E., Clamp, M. & Durbin, R. GeneWise and genomewise. Genome Res. 14, 988–995 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Edgar, R. C. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
McKenna, A. et al. The genome analysis toolkit: a mapreduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156–2158 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).Article 

Google Scholar 
Robinson, J. D., Bunnefeld, L., Hearn, J., Stone, G. N. & Hickerson, M. J. ABC inference of multi-population divergence with admixture from unphased population genomic data. Mol. Ecol. 23, 4458–4471 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Nazareno, A. G., Bemmels, J. B., Dick, C. W. & Lohmann, L. G. Minimum sample sizes for population genomics: an empirical study from an Amazonian plant species. Mol. Ecol. Resour. 17, 1136–1147 (2017).CAS 
PubMed 
Article 

Google Scholar 
Willing, E. M., Dreyer, C. & van Oosterhout, C. Estimates of genetic differentiation measured by FST do not necessary require large sample size when using many SNP markers. PLoS One 7, e2649 (2012).Article 
CAS 

Google Scholar 
Keenan, K., Mcginnity, P., Cross, T. F., Crozier, W. W. & Prodöhl, P. A. diveRsity: an Rpackage for the estimation and exploration of population genetics parameters and their associated errors. Methods Ecol. Evol. 4, 782–788 (2013).Article 

Google Scholar 
Rellstab, C., Gugerli, F., Eckert, I. A., Hancock, M. A. & Holderegger, R. A practical guide to environmental assocaition analysis in landscape genomics. Mol. Ecol. 24, 4348–4370 (2015).PubMed 
Article 

Google Scholar 
Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly 6, 80–92 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Xie, C. et al. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316–W322 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Araújo, M. B. & New, M. Ensemble forecasting of species distributions. Trends Ecol. Evol. 22, 42–47 (2007).PubMed 
Article 

Google Scholar 
Araújo, M. B. et al. Standards for distribution models in biodiversity assessments. Sci. Adv. 5, eaat4858 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Boria, R. A., Olson, L. E., Goodman, S. M. & Anderson, R. P. Spatial filtering to reduce sampling bias can improve the performance of ecological niche models. Ecol. Model. 275, 73–77 (2014).Article 

Google Scholar 
Anderson, R. P. & Raza, A. The effect of the extent of the study region on GISmodels of species geographic distributions and estimates of niche evolution: preliminary tests with montane rodents (genus Nephelomys) in Venezuela. J. Biogeogr. 37, 1378–1393 (2010).Article 

Google Scholar 
Pearson, R. G., Raxworthy, C., Nakamura, M. & Peterson, A. T. Predicting species distributions from small numbers of occurrence records: a test case using crypticgeckos in Madagascar. J. Biogeogr. 34, 102–117 (2007).Article 

Google Scholar 
Heming, N. M., Dambros, C. & Gutiérrez, E. E. ENMwizard: advanced techniques for Ecological Niche Modeling made easy. https://github.com/HemingNM/ENMwizard (2018).Naimi, B., Hamm, N. A. S., Groen, T. A., Skidmore, A. K. & Toxopeus, A. G. Where is positional uncertainty a problem for species distribution modelling. Ecography 37, 191–203 (2014).Article 

Google Scholar 
Muscarella, R. et al. ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for MAXENT ecological niche models. Methods Ecol. Evol. 5, 1198–1205 (2014).Article 

Google Scholar 
Owens, H. L. et al. Constraints on interpretation of ecological niche models by limited environmental ranges on calibration areas. Ecol. Model. 263, 10–18 (2013).Article 

Google Scholar 
Akaike, H. New look at statistical-model identification. IEEE Trans. Autom. Control AC19, 716–723 (1974).ADS 
MathSciNet 
MATH 
Article 

Google Scholar 
Phillips, S. J. et al. Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecol. Appl. 19, 181–197 (2009).PubMed 
Article 

Google Scholar 
Barbet-Massin, M., Jiguet, F., Albert, C. H. & Thuiller, W. Selecting pseudo-absences for species distribution models: how, where and how many? Methods Ecol. Evol. 3, 327–338 (2012).Article 

Google Scholar 
Bellard, C. et al. Will climate change promote future invasions? Glob. Change Biol. 19, 3740–3748 (2013).ADS 
Article 

Google Scholar 
Elith, J., Kearney, M. & Phillips, S. The art of modelling range-shifting species. Methods Ecol. Evol. 1, 330–342 (2010).Article 

Google Scholar 
Anantharaman, R., Hall, K., Shah, V. B. & Edelman, A. Circuitscape in Julia: high performance connectivity modelling to support conservation decisions. Proc. JuliaCon Conf. 1, 58 (2020).Article 

Google Scholar 
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar 
Anderson, D. R. & Burnham, K. P. Avoiding pitfalls when using information-theoretic methods. J. Wildl. Manag. 66, 912–918 (2002).Article 

Google Scholar 
Van Strien, M. J., Keller, D. & Holderegger, R. A new analytical approach to landscape genetic modelling: least-cost transect analysis and linear mixed models. Mol. Ecol. 21, 4010–4023 (2012).Article 

Google Scholar 
Bartoń, K. MuMIn: multi-model inference, R package version 1.9.13 (2013).Zhang, G. et al. Comparative genomics reveal insights into avian genome evolution and adaptation. Science 346, 1311–1320 (2014).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Roesti, M., Kueng, B., Moser, D. & Berner, D. The genomics of ecological vicariance in threespine stickleback fish. Nat. Commun. 6, 8767 (2015).ADS 
CAS 
PubMed 
Article 

Google Scholar  More