More stories

  • in

    Phenotypic trait variation in a long-term multisite common garden experiment of Scots pine in Scotland

    Seed sampling and germinationSeed from ten trees from each of 21 native Scottish Scots pine populations (Table 1) were collected in March 2007 and germinated at the James Hutton Institute, Aberdeen (latitude 57.133214, longitude −2.158764) in June 2007. Populations were chosen to represent the species’ native range in Scotland and to include three populations from each of the seven seed zones (Fig. 2). There was no selection of seed-trees on the basis of any traits except for the possession of cones on the date of sampling. Ten seed trees were sampled from each population according to a spatial protocol designed to cover a circle of approximately 1 km in diameter located around a central tree. The sampling strategy identified nine points each in a pre-determined random direction from the central point, whilst stratifying the number sampled with increasing distance from the central point in the ratio 1: 3: 5. This strategy avoids over-sampling the areas close to the centre point. For smaller fragments of woodland, or where only a few trees with cones were present, then the directions of the sampled trees from the central tree were maintained to give a wide coverage of the woodland area, but the distances between trees varied but were never closer than 50 m. To break dormancy, seeds were soaked for 24 hours on the benchtop at room temperature, after which they were stored in wet paper towels and refrigerated in darkness at 3–5 °C for approximately 4 weeks. Seeds were kept moist and transferred to room temperature until germination began (approx. 5–7 days), then transplanted to 8 cm × 8 cm × 9 cm, 0.4 L pots filled with Levington’s C2a compost and 1.5 g of Osmocote Exact 16–18 months slow release fertiliser and kept in an unheated glasshouse. The compost was covered with a layer of grit to reduce moss and liverwort growth. Seedlings from the same mother tree are described as a family and are assumed to be half-siblings.Table 1 Locations and basic environmental data for the populations sampled for seed to establish the trial. See the maternal traits dataset15 for precise data for each mother tree sampled.Full size tableExperimental design: nurseriesThe full collection consisted of 210 families (10 families from each of 21 populations) each consisting of 24 half sibling progeny (total 5,040 individuals); needle tissue was sampled from each seedling and preserved for long term storage, one needle on silica gel, 2–5 needles at −20 °C. After transfer into pots, 8 seedlings per family were moved to one of three nurseries (total 1,680 seedlings per nursery): outdoors at Inverewe Gardens in western Scotland (nursery in the west of Scotland: coded NW, latitude 57.775714, longitude −5.597181, Fig. 2); outdoors in a fruit cage (to minimise browsing) at the James Hutton Institute in Aberdeen (nursery in the east of Scotland: NE); in an unheated glasshouse at the James Hutton Institute in Aberdeen (nursery in a glasshouse: NG). Trees were arranged in 40 randomised trays (blocks) in each nursery. Each block contained two trees per population (total 42 trees). Watering was automatic in NG, and manually as required for NE and NW. No artificial light was used in any of the nurseries. In May 2010 the seedlings from NG were moved outdoors to Glensaugh in Aberdeenshire (latitude 56.893567, longitude −2.535736). In 2010 all plants were repotted into 19 cm (3 L) pots containing Levingtons CNSE Ericaceous compost with added Osmocote STD 16–18 month slow release fertilizer.Experimental design: field sitesIn 2012 the trees were transplanted to one of three field sites: Yair in the Scottish Borders (field site in the south of Scotland: FS, latitude 55.603625, longitude −2.893025); Glensaugh (field site in the east of Scotland: FE); and Inverewe (field site in the west of Scotland: FW). All trees transplanted to FS were raised in the NG and all but four of the trees transplanted to FE were raised locally in the NE (the remainder were grown in NG). In contrast, following mortality and ‘beating up’ (filling gaps where saplings had died), the FW experiment ultimately contained cohorts of trees raised in each of the three nurseries as follows: 290 grown locally in the NW; 132 were grown in the NG; and 82 were grown in the NE.Site historiesThe Yair site (FS) had previously been used for growing Noble fir (Abies procera) for Christmas trees and Lodgepole pine (Pinus contorta), a section of the former were felled and chipped to create a clear area prior to planting. The planting site is also adjacent to a large block of commercial Sitka spruce (Picea sitchensis) forestry, and the Glenkinnon Burn Site of Special Scientific Interest (SSSI NatureScot site code 736; EU site code 135445), an area of mixed broadleaf woodland. Prior to planting, major areas of tall weeds were strimmed. The site was protected by a deer fence. The experiment was planted 8–11 October 2012. The Glensaugh site (FE) is in Forestry Compartment 3 of the Glensaugh Research Station, adjacent to Cleek Loch. It is thought to have been cleared of Scots pine and Larch (Larix decidua) around 1917, after which it reverted to rough grazing. An attempt to reseed part of the site in the 1980s was unsuccessful and it quickly reverted to rough grazing for a second time. The whole site (within which the experimental area is embedded) was deer fenced and re-planted under the Scottish Rural Development Programme (SRDP) in 2012. The experimental plot was planted up 7–9 March 2012. The Inverewe site (FW) had previously been a Sitka spruce and Lodgepole pine plantation (50:50 mix) that had been clear-felled in 2010 following substantial windthrow. The experimental site was deer fenced in early 2012, and the experiment was planted 12–16 March 2012, followed by beating up on 27–28 March 2013 and 22–24 October 2013. There had been minimal preparation of the site in line with current practice for restocking sites. The experimental site is included in the Inverewe Forest Plan, which included deer fencing of a larger area (2014) around the experimental site. Planting of this area was completed in early 2015, funded by NTS (National Trust for Scotland), although natural regeneration is also taking place.At each site, trees were planted in randomised blocks at 3 m × 3 m spacing. There are four randomised blocks in both FS and FE and three in FW. A guard row of Scots pine trees was planted around the periphery of the blocks and between blocks B and C at FS. Each block comprised one individual from each of eight (of the 10 sampled) families per 21 populations (168 trees). Although most families (N = 159) were represented at each of the three sites, families with insufficient trees (N = 9) were replaced in one site (FS) with a different family from the same population. Each experimental site was designed with redundancy such that, if thinning becomes necessary as the trees mature, then the systematic removal of trees (i.e. trees 1,3,5,7, etc of row 1, and 2,4,6,8, etc of row 2, 1,3,5,7,etc of row 3) will maintain a balanced design of the experiment, with sufficient family and population representation to provide an ongoing experiment with full geographic coverage.The field sites generally experience different climates, with FW typically warmer and wetter and with more growing degree days per year and a much longer growing season than both FE and FS (Table 2). The coldest site with the shortest growing season is generally FE.Table 2 Average climatic variables at field sites Glensaugh (FE), Inverewe (FW) and Yair (FS) from planting in 2012 until 2020. Climatic variables are derived from data provided by the Met Office (daily mean, minimum and maximum temperatures and monthly rainfall).Full size tablePhenotype assessmentsMaternal traitsFollowing seed collection, a range of traits were measured in the mother trees in order to control for maternal effects in subsequent measurements of their progeny (Table 3). For each mother tree, measurements of height and diameter at breast height (DBH) were taken, and ten cones were collected and assessed in detail. Cone width and length were measured prior to drying the cones (when they were still closed). Cone weight was measured post-drying. Seed removed from each cone was assessed for total weight (after wings had been removed) and for the count and percentage of seeds which were classed as viable (viable seed were those which had both a wing and an obvious seed). No further seed sorting was applied.Table 3 Traits assessed in mother trees, cones, seeds (dataset: Maternal), nursery seedlings (dataset: Nursery) and field trials (dataset: Field). Within the datasets, traits are recorded in a single column for each year using the format Code-Year (e.g. absolute height in 2008 = HA08) except for the maternal traits datasets which were all measured in 2007.Full size tableNursery traitsSeedling phenotype assessments were performed annually from 2007–2010 for three different trait types: phenology (budburst and growth cessation); form (total number of buds, needle length); cumulative growth (stem diameter and height, canopy width). Measurements of tree form and cumulative growth traits were taken after the end of each growing season. Phenology was assessed weekly during the spring and autumn of 2008 for budburst and growth cessation, respectively. Budburst was defined as the number of days from 31 March 2008 to the time when newly emerged green needles were observed (budburst stage 6: Fig. 3). In some seedlings in 2008, a secondary flush of growth occurred from terminal buds that had formed during the summer of that year. Therefore, growth cessation was defined retrospectively as the number of days from 10 September 2008 to the date when a terminal bud had formed on the leading shoot of the seedling, providing no further growth was observed either on the stem below that bud, nor from the bud itself. Canopy width (widest point) was measured at two perpendicular points in the horizontal plane. Needle length was measured for three needles per tree. Mortality was recorded each year from 2007 to 2010.Fig. 3Phenological stages of bud burst in Pinus sylvestris assessed in field trials. Inset numbers correspond to budburst stage. Budburst stage 1: bud dormant; 2: bud swelling and showing signs of linear expansion; 3: scales open at base revealing green tissue. Remaining bud remains encased by smooth bud scales; 4: scales open along length of shoot revealing green tissue and partially visible needles; 5: white tipped needles visible along length of the shoot; 6: green needles emerging away from the shoot (bottle brush appearance) along its entire length; 7: Needles have separated and next year’s terminal bud is usually formed and clearly visible.Full size imageField traitsTree height was measured in the field in the winter after each growing season from 2013 at FE and FW, and from 2014 to 2020 at all sites. Height was taken as the vertical measurement in cm from top bud straight to the ground. Basal stem diameter was measured at the end of the growing season for trees growing at FE and FW from 2014 to 2020 and for FS in 2020.Phenology assessments were performed in spring at each site from 2015 to 2019. Seven distinct stages of budburst (assessed on the terminal bud) were defined (Fig. 3) although only stages 4 to 6 are included in the dataset and considered for analysis due to high proportions of missing data for the early and late stages. Each tree was assessed for budburst stage at weekly intervals from early spring until budburst was complete. In order to allow comparisons within and among sites and years, the date at which each stage of budburst occurred was considered relative to 31 March of that year. For example, 25 May 2019 is recorded as 55 days since 31 March 2019. The duration of budburst (time taken to reach stage 6 from stage 4) was also estimated.When trees progressed through budburst stages rapidly, skipping a stage between assessments, a mean value was taken from the two assessment dates. For example, if a tree was at stage 4 on day 55 and was recorded as stage 6 at the next assessment on day 62, it is assumed to have reached stage 5 at day 58.5. More

  • in

    Longitudinal analysis of the Five Sisters hot springs in Yellowstone National Park reveals a dynamic thermoalkaline environment

    Mueller, R. C. et al. An emerging view of the diversity, ecology, and function of Archaea in alkaline hydrothermal environments. FEMS Microbiol. Ecol. 97, fiaa246 (2020).
    Google Scholar 
    López-López, O., Cerdán, M.-E. & González-Siso, M.-I. Thermus thermophilus as a source of thermostable lipolytic enzymes. Microorganisms 3, 792–808 (2015).PubMed 
    PubMed Central 

    Google Scholar 
    Sahay, H. et al. Hot springs of Indian Himalayas: Potential sources of microbial diversity and thermostable hydrolytic enzymes. 3 Biotech 7, 118 (2017).PubMed 
    PubMed Central 

    Google Scholar 
    Patel, A. K., Singhania, R. R., Sim, S. J. & Pandey, A. Thermostable cellulases: Current status and perspectives. Bioresour Technol 279, 385–392 (2019).CAS 
    PubMed 

    Google Scholar 
    Decastro, M.-E., Rodríguez-Belmonte, E. & González-Siso, M.-I. Metagenomics of thermophiles with a focus on discovery of novel thermozymes. Front. Microbiol. 7, 1521–1521 (2016).PubMed 
    PubMed Central 

    Google Scholar 
    Meslé, M. M. et al. Isolation and characterization of lignocellulose-degrading geobacillus thermoleovorans from Yellowstone National Park. Appl. Environ. Microbiol. 88, e0095821 (2022).PubMed 

    Google Scholar 
    Verma, P., Yadav, A. N., Shukla, L., Saxena, A. K. & Suman, A. Hydrolytic enzymes production by thermotolerant Bacillus altitudinis IARI-MB-9 and Gulbenkiania mobilis IARI-MB-18 isolated from Manikaran hot springs. Int. J. Adv. Res. 3, 1241–1250 (2015).CAS 

    Google Scholar 
    Wu, B. et al. Microbial sulfur metabolism and environmental implications. Sci. Total Environ. 778, 146085 (2021).ADS 
    CAS 
    PubMed 

    Google Scholar 
    Lavrentyeva, E. V. et al. Bacterial diversity and functional activity of microbial communities in hot springs of the Baikal Rift Zone. Microbiology 87, 272–281 (2018).CAS 

    Google Scholar 
    Miller Scott, R., Strong Aaron, L., Jones Kenneth, L. & Ungerer Mark, C. Bar-Coded pyrosequencing reveals shared bacterial community properties along the temperature gradients of two alkaline hot springs in Yellowstone National Park. Appl. Environ. Microbiol. 75, 4565–4572 (2009).ADS 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Sharp, C. E. et al. Humboldt’s spa: Microbial diversity is controlled by temperature in geothermal environments. ISME J. 8, 1166–1174 (2014).CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Stefanova, K. et al. Archaeal and bacterial diversity in two hot springs from geothermal regions in Bulgaria as demostrated by 16S rRNA and GH-57 genes. Int. Microbiol. 18, 217–223 (2015).CAS 
    PubMed 

    Google Scholar 
    Hou, W. et al. A comprehensive census of microbial diversity in hot springs of Tengchong, Yunnan Province China using 16S rRNA gene pyrosequencing. PLoS ONE 8, e53350 (2013).ADS 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Sahm, K. et al. High abundance of heterotrophic prokaryotes in hydrothermal springs of the Azores as revealed by a network of 16S rRNA gene-based methods. Extremophiles 17, 649–662 (2013).CAS 
    PubMed 

    Google Scholar 
    Purcell, D. et al. The effects of temperature, pH and sulphide on the community structure of hyperthermophilic streamers in hot springs of northern Thailand. FEMS Microbiol. Ecol. 60, 456–466 (2007).CAS 
    PubMed 

    Google Scholar 
    Meyer-Dombard, D. R. & Amend, J. P. Geochemistry and microbial ecology in alkaline hot springs of Ambitle Island, Papua New Guinea. Extremophiles 18, 763–778 (2014).CAS 
    PubMed 

    Google Scholar 
    de Leon, K. B., Gerlach, R., Peyton, B. M. & Fields, M. W. Archaeal and bacterial communities in three alkaline hot springs in Heart Lake Geyser Basin, Yellowstone National Park. Front. Microbiol. 4, 10 (2013).
    Google Scholar 
    Boomer, S. M., Noll, K. L., Geesey, G. G. & Dutton, B. E. Formation of multilayered photosynthetic biofilms in an alkaline thermal spring in Yellowstone National Park, Wyoming. Appl. Environ. Microbiol. 75, 2464–2475 (2009).ADS 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Wang, S. et al. Greater temporal changes of sediment microbial community than its waterborne counterpart in Tengchong hot springs, Yunnan Province, China. Sci. Rep. 4, 7479 (2014).ADS 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Sun, Y., Liu, Y., Pan, J., Wang, F. & Li, M. Perspectives on cultivation strategies of archaea. Microb. Ecol. 79, 770–784 (2020).PubMed 

    Google Scholar 
    Brock, T. D. Life at high temperatures. Science 158, 1012 (1967).ADS 
    CAS 
    PubMed 

    Google Scholar 
    Christiansen, R. L. The Quaternary and Pliocene Yellowstone Plateau volcanic field of Wyoming, Idaho, and Montana. Professional Paper (2001).Rowe, J. J., Fournier, R. & Morey, G. Chemical analysis of thermal waters in Yellowstone National Park, Wyoming, 1960–65. USGS https://doi.org/10.3133/b1303 (1973).Article 

    Google Scholar 
    Fournier, R., Thompson, M. J. & Hutchinson, R. A. The geochemistry of hot spring waters at Norris Geyser Basin, Yellowstone National Park. International symposium on water-rock interactions (1992).Podar, P. T., Yang, Z., Björnsdóttir, S. H. & Podar, M. Comparative analysis of microbial diversity across temperature gradients in hot springs from Yellowstone and Iceland. Front. Microbiol. 11, 1625 (2020).PubMed 
    PubMed Central 

    Google Scholar 
    Pala, C. et al. Environmental drivers controlling bacterial and archaeal abundance in the sediments of a Mediterranean lagoon ecosystem. Curr. Microbiol. 75, 1147–1155 (2018).CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Foyer, C. H., Noctor, G. & Hodges, M. Respiration and nitrogen assimilation: Targeting mitochondria-associated metabolism as a means to enhance nitrogen use efficiency. J. Exp. Bot. 62, 1467–1482 (2011).CAS 
    PubMed 

    Google Scholar 
    Ershanovich, V. N. et al. Nitrogen assimilation enzymes in Bacillus subtilis mutants with hyperproduction of riboflavin. Mol. Gen. Mikrobiol. Virusol. 2005(3), 29–34 (2005).
    Google Scholar 
    Offre, P., Spang, A. & Schleper, C. Archaea in biogeochemical cycles. Annu Rev Microbiol 67, 437–457 (2013).CAS 
    PubMed 

    Google Scholar 
    Cabello, P., Roldán, M. D. & Moreno-Vivián, C. Nitrate reduction and the nitrogen cycle in archaea. Microbiology 150, 3527–3546 (2004).CAS 
    PubMed 

    Google Scholar 
    Graupner, M., Xu, H. & White, R. H. The pyrimidine nucleotide reductase step in riboflavin and F(420) biosynthesis in archaea proceeds by the eukaryotic route to riboflavin. J. Bacteriol. 184, 1952–1957 (2002).CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Chernyh, N. A. et al. Dissimilatory sulfate reduction in the archaeon “Candidatus Vulcanisaeta moutnovskia” sheds light on the evolution of sulfur metabolism. Nat. Microbiol. 5, 1428–1438 (2020).CAS 
    PubMed 

    Google Scholar 
    Castelle, C. J. & Banfield, J. F. Major new microbial groups expand diversity and alter our understanding of the tree of life. Cell 172, 1181–1197 (2018).CAS 
    PubMed 

    Google Scholar 
    Williams, T. A. et al. Integrative modeling of gene and genome evolution roots the archaeal tree of life. Proc. Natl. Acad. Sci. U.S.A. 114, E4602–E4611 (2017).ADS 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Guy, L. & Ettema, T. J. G. The archaeal ‘TACK’ superphylum and the origin of eukaryotes. Trends Microbiol. 19, 580–587 (2011).CAS 
    PubMed 

    Google Scholar 
    Wang, Y., Wegener, G., Hou, J., Wang, F. & Xiao, X. Expanding anaerobic alkane metabolism in the domain of Archaea. Nat. Microbiol. 4, 595–602 (2019).CAS 
    PubMed 

    Google Scholar 
    Hedlund, B. P. et al. Uncultivated thermophiles: Current status and spotlight on ‘Aigarchaeota’. Curr. Opin. Microbiol. 25, 136–145 (2015).CAS 
    PubMed 

    Google Scholar 
    Reichart, N. J. et al. Activity-based cell sorting reveals responses of uncultured archaea and bacteria to substrate amendment. ISME J. 14, 2851–2861 (2020).CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Hua, Z.-S. et al. Genomic inference of the metabolism and evolution of the archaeal phylum Aigarchaeota. Nat. Commun. 9, 2832 (2018).ADS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Beam, J. P. et al. Ecophysiology of an uncultivated lineage of Aigarchaeota from an oxic, hot spring filamentous “streamer” community. ISME J. 10, 210–224 (2016).CAS 
    PubMed 

    Google Scholar 
    Gonsior, M. et al. Yellowstone hot springs are organic chemodiversity hot spots. Sci. Rep. 8, 14155 (2018).ADS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Gibson, M. L. & Hinman, N. W. Mixing of hydrothermal water and groundwater near hot springs, Yellowstone National Park (USA): Hydrology and geochemistry. Hydrogeol. J. 21, 919–933 (2013).ADS 
    CAS 

    Google Scholar 
    Campbell, K. M. et al. Sulfolobus islandicus meta-populations in Yellowstone National Park hot springs. Environ. Microbiol. 19, 2334–2347 (2017).PubMed 

    Google Scholar 
    Thiel, V. et al. The dark side of the mushroom spring microbial mat: Life in the shadow of chlorophototrophs. I. Microbial diversity based on 16S rRNA gene amplicons and metagenomic sequencing. Front. Microbiol. 7, 919 (2016).PubMed 
    PubMed Central 

    Google Scholar 
    Caporaso, J. G. et al. Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc. Natl. Acad. Sci. U.S.A. 108, 4516–4522 (2011).ADS 
    CAS 
    PubMed 

    Google Scholar 
    Parada, A. E., Needham, D. M. & Fuhrman, J. A. Every base matters: Assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ. Microbiol. 18, 1403–1414 (2016).CAS 
    PubMed 

    Google Scholar 
    Apprill, A., McNally, S., Parsons, R. & Weber, L. Minor revision to V4 region SSU rRNA 806R gene primer greatly increases detection of SAR11 bacterioplankton. Aquat. Microb. Ecol. 75, 129–137 (2015).
    Google Scholar 
    Thompson, L. R. et al. A communal catalogue reveals Earth’s multiscale microbial diversity. Nature 555, 457–463 (2017).ADS 

    Google Scholar 
    Eloe-Fadrosh, E. A., Ivanova, N. N., Woyke, T. & Kyrpides, N. C. Metagenomics uncovers gaps in amplicon-based detection of microbial diversity. Nat. Microbiol. 1, 15032 (2016).CAS 
    PubMed 

    Google Scholar 
    Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26, 2460–2461 (2010).CAS 
    PubMed 

    Google Scholar 
    Edgar, R. C. UNOISE2: Improved error-correction for Illumina 16S and ITS amplicon sequencing. BioRxiv https://doi.org/10.1101/081257 (2016).Article 

    Google Scholar 
    Murali, A., Bhargava, A. & Wright, E. S. IDTAXA: A novel approach for accurate taxonomic classification of microbiome sequences. Microbiome 6, 140 (2018).PubMed 
    PubMed Central 

    Google Scholar 
    Parks, D. H. et al. A complete domain-to-species taxonomy for Bacteria and Archaea. Nat. Biotechnol. 38, 1079–1086 (2020).CAS 
    PubMed 

    Google Scholar 
    Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Mol. Biol. Evol. 30, 772–780 (2013).CAS 
    PubMed 
    PubMed Central 

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

    Google Scholar 
    Letunic, I. & Bork, P. Interactive Tree Of Life (iTOL) v5: An online tool for phylogenetic tree display and annotation. Nucleic Acids Res. 49, W293–W296 (2021).CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Matsen, F. A., Kodner, R. B. & Armbrust, E. V. pplacer: Linear time maximum-likelihood and Bayesian phylogenetic placement of sequences onto a fixed reference tree. BMC Bioinform. 11, 538 (2010).
    Google Scholar 
    Chong, J., Liu, P., Zhou, G. & Xia, J. Using MicrobiomeAnalyst for comprehensive statistical, functional, and meta-analysis of microbiome data. Nat. Protoc. 15, 799–821 (2020).CAS 
    PubMed 

    Google Scholar 
    Chambers, M. C. et al. A cross-platform toolkit for mass spectrometry and proteomics. Nat. Biotechnol. 30, 918–920 (2012).CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Pluskal, T., Castillo, S., Villar-Briones, A. & Oresic, M. MZmine 2: Modular framework for processing, visualizing, and analyzing mass spectrometry-based molecular profile data. BMC Bioinform. 11, 395 (2010).
    Google Scholar 
    Patiny, L. & Borel, A. ChemCalc: A building block for tomorrow’s chemical infrastructure. J. Chem. Inf. Model. 53, 1223–1228 (2013).CAS 
    PubMed 

    Google Scholar 
    Chong, J., Wishart, D. S. & Xia, J. Using MetaboAnalyst 4.0 for comprehensive and integrative metabolomics data analysis. Curr. Protoc. Bioinform. 68, e86 (2019).
    Google Scholar 
    Liu, G., Lee, D. P., Schmidt, E. & Prasad, G. L. Pathway analysis of global metabolomic profiles identified enrichment of caffeine, energy, and arginine metabolism in smokers but not moist snuff consumers. Bioinform. Biol. Insights 13, 1177932219882961–1177932219882961 (2019).PubMed 
    PubMed Central 

    Google Scholar 
    Xia, J. & Wishart, D. S. MetPA: A web-based metabolomics tool for pathway analysis and visualization. Bioinformatics 26, 2342–2344 (2010).CAS 
    PubMed 

    Google Scholar 
    Huber, W. et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat. Methods 12, 115–121 (2015).CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Rohart, F., Gautier, B., Singh, A. & Lé Cao, K.-A. mixOmics: An R package for ’omics feature selection and multiple data integration. PLoS Comput. Biol. 13, e1005752–e1005752 (2017).ADS 
    PubMed 
    PubMed Central 

    Google Scholar  More

  • in

    The dominant mesopredator and savanna formations shape the distribution of the rare northern tiger cat (Leopardus tigrinus) in the Amazon

    Most records of N-tiger cats were from savanna environments, and it was not surprising that this vegetative formation has a key influence on the N-tiger cat range in the Amazon. The bulk of the L. t. tigrinus distribution lies in the savannas, dry forests and shrublands of the Cerrado and Caatinga biomes. These are also the areas with the vast majority of records for this lowland subspecies (Supplemental Fig. S5). Hence, L. t. tigrinus is more associated with savannas and savanna-like environments than with rainforests. In fact, more than 80% of the records in the Amazon were within 100 km of a savanna patch. Colonization of the northern savanna formations of the Amazon by the N-tiger cat likely occurred during the forest-savanna shifts of the glacial period18, and the cat currently shows a patchy distribution. Strong evidence of established biogeographic corridor connections between the savannas of the Cerrado and those of the Amazon exists, suggesting northward expansion of the former during glacial periods, perhaps predating the Last Glacial Maximum19,20,21. Further corroborating this evidence, tiger cat ‘gene flow’ niche modelling showed prior connectivity between the Guiana population and that of Central Brazil and no connectivity with the Andean population22. Additionally, Guianan tiger cat skin patterns are found in savanna and transitional savanna/Amazon areas and in the semiarid shrub-woodland of Brazil and are very distinct from the patterns of the tiger cats from the Andes of northwestern South America and Central America (Supplementary Information Fig. S6).The bioclimatic variables in the best model also supported the cat’s preference for savanna areas. The best model indicated a positive effect of precipitation in the driest month on the probability of the presence of the N-tiger cat, likely indicating the Aw/As climates of tropical savannas23. These climates are marked by seasonal variation in rainfall, with a pronounced dry season. Higher rainfall during the dry season favors the growth of vegetation, which results in some tree cover within the savannas. Thus, our results agree with previous research suggesting that tiger cats avoid open savanna formations24. Similarly, the species had a significant negative response to net primary productivity. This also supports the species’ avoidance of dense lowland rainforests, which are the most productive habitats. In the Amazon biome, the least productive areas are found in more open landscapes25.The N-tiger cat’s range considered from an ecoregion perspective12 could biogeographically explain its distribution in the Amazon. All records but 2 fell within Guiana savannas, Guiana highland forest, Guiana rainforest, part of the Uatumã-Trombetas rainforest bordering the Guianas or all of it connecting to Gurupá and Monte Alegre varzea forests, as well as Marajó varzeas, the interfluve Tocantins-Araguaia/Maranhão, and the southern block of the interfluve Xingu/Tocantins-Araguaia. There were two records from the Negro-Branco moist forest, which also includes savanna-like “campinarana” formations. The range also reaches the transitional babaçu palm forests of Maranhão and the Mato Grosso seasonal forests (Supplementary Information Fig. S7, Table S3). The N-tiger cat’s range in the Amazon was determined by combining records with species distribution modeling, also matching the ecoregion perspective.Outside the Guiana Shield and likely the savanna patches of the region of the Upper Negro River, in other parts of the Amazon, the N-tiger cat seems to be restricted to the forests of the eastern Amazon, along the arc of deforestation and to transitional areas with savanna formations. The presence and absence points at camera-trapping sites could explain the N-tiger cat’s range in the Amazon and define its distribution range in the biome. Absence points, for instance, were usually located in dense rainforest habitats throughout the Amazon biome.The species may occasionally occupy rainforests, such as those of the Guianas, where it tends to be very rare. At a site in central Suriname, after an enormous trapping effort of  > 20,000 trap days in four years by cat specialists, over an area  > 1100 km2, no records of the N-tiger cat were found (Supplementary Information Table S2), although its presence is expected in that area26. This finding attests to the inherent rarity of this felid in its limited range within the Amazon. However, could its association with the arc of deforestation be related to the replacement of forest by bushy savanna-like vegetation that succeeds abandoned pastures? The other currently recognized subspecies, L. t. pardinoides (the Andean tiger cat) and L. t. oncilla (the oncilla), and the recently split southern tiger cat L. guttulus are all associated with forested areas. Conversely, L. t. tigrinus has higher abundance and is mostly found in the nonforested habitats of the Cerrado and Caatinga domains of Brazil and only rarely in rainforests. Thus, L. t. tigrinus may be an open-habitat (sub)species. However, within savannas, N-tiger cats are restricted to denser savanna formations, with open savannas deemed unsuitable24. In the semiarid Caatinga, the N-tiger cat also prefers denser formations27,28.One of the most interesting findings was the clear relationship between the ranges of the dominant mesopredator and subordinate species. The ranges of ocelots and N-tiger cats in the Amazon were diametrically opposite (Fig. 1), a finding never recorded for felids. The reported ocelot densities and relative abundance indexes (RAIs) in the Amazon range from 0.29 to 0.95 ind/km2 and 0.07–13.2 ind/100 trap-days, respectively7,29. Thus, the expected ocelot density found using modeling that allows for N-tiger cat presence is very low (Fig. 2A). In the Rupununi, the ocelot:N-tiger cat RAI ratio was roughly 10:1, with a very low RAI and expected density for N-tiger cats (see Supplementary Material). The only other relative abundance estimate of tiger cats presented for the Amazon30 was not confirmed as an estimate of tiger cats following inspection of the original records by the authors but as an estimate of margays or ocelots. This antagonistic relationship between ocelots and all other small cat species in their area of sympatry is quite impressive. It is density-dependent, as it seems to take effect only above an ocelot density threshold of 0.12 ind./km231. The influence can range from patterns of density, distribution, and occupancy to spatial and temporal use. Conversely, such an impact was not detected when either the small cats or ocelots were compared to the larger cats31,32,33,34,35.In view of the Red List assessments and applying the limited estimates presented, the expected total population size for N-tiger cats in the Amazon would be approximately 150 and 1622 individuals, considering their AOO or EOO, respectively. Applying the IUCN’s formula for mature individuals8, these numbers would be 45 and 487 individuals for the AOO and EOO, respectively.The ocelot’s preference for very dense rainforests may explain the low probability of N-tiger cat occurrence within the Amazon biome. Notably, most tiger cat records from rainforests and all those from premontane forests came from the Guiana Shield, a region where tropical grasslands and savannas dot more forested landscapes. The Guiana Highlands and Pantepui ecoregions, which make up a considerable portion of the shield, tend to have low ocelot densities (below 0.30 ind/km2), although they do contain some rainforest. Ocelot densities reach some of their lowest values in the Guianan savanna ecoregion (mean ocelot density of 0.029 in the savanna formations), where the N-tiger cat probability of occurrence was highest. At the Karanambu site in the Rupununi, all ocelot records came from either gallery forests or forest patches embedded in the savanna. Although the data did not allow us to test further hypotheses, it is likely that spatial partitioning occurs in the Guiana Shield, with N-tiger cats favoring habitats that are more open. Conversely, areas farther west in the Amazon biome, other than the predicted area, do not have any major savanna patches and are covered mostly by lowland tropical rainforest formations, where ocelots can potentially reach densities in excess of 0.7 ind/km2. Of all Amazonian records of N-tiger cats, only one came from west of the 68th meridian: a preserved specimen from Puerto Leguizamo on the Putumayo River in Colombia. The specimen was identified as L. t. pardinoides by its collector, so it most likely represents an individual that came down from the foothills of the Andes. Alternatively, it could have been caught in the Andean foothills but labeled generally as from Puerto Leguizamo, as museum records do not always present precise locations, like most of those from our dataset; thus, they could represent a broader region, not a single collection location.The records of L. t. tigrinus in the Monte-Alegre Várzea ecoregion and Tapajós-Xingu Moist Forest ecoregion (which shares a border with the Amazon River) are actually from the small savanna patches of Terra Santa and Alter do Chão, respectively, which are imbedded within the forests of these ecoregions. Similarly, the Negro-Branco Moist Forest ecoregion includes open-canopy white sand forests with savanna-like vegetation, known as ‘campinaranas’36.Although our model predicted a high probability of N-tiger cat presence in the Marajó Várzea ecoregion, the records from the island came from savanna patches and not from flooded forests and mangroves. Hence, we did not include such large areas in the AOO for the subspecies. It is likely that the highly predicted probability of presence there is an artifact of low predicted ocelot density. Nevertheless, the environment there is not suitable for either cat. Our ocelot density model was highly significant and explained almost 50% of the variation in ocelot density. The remaining variation was related to either other variables that could not be measured via satellite imagery (such as prey availability) or the sampling design of the different studies. Nonetheless, ocelot densities predicted from our model across the Amazon were within the expected range for the species29.Why are N-tiger cats absent in camera-trapping studies in Amazonian forests throughout the biome? The most straightforward answer seems to be because they simply are not there (central and western Amazon) or, where present, their numbers are extremely low (Guianas and eastern Amazon). The lack of surveys cannot be cited as a potential reason for their apparent absence because the studies that did not detect the species were conducted throughout the Amazon biome, in all nine Amazonian countries. Some of the areas have been surveyed for several years—or decades in some cases—and have failed to record a single individual (Supplementary Information Table S2). Typically, N-tiger cats appear, even prominently, on cameras in other biomes, such as in the savannas of the Cerrado and semiarid scrub of the Caatinga domain in Brazil, including sites where ocelots are present24,27,37. Clouded tiger cats (L. t. pardinoides) have also been frequently recorded on cameras in the Andes, higher than 1500 m above sea level34,38, but not in lowland Amazonian forests. This finding indicates that the N-tiger cat is not camera-shy. In northern Brazilian savannas, its density can reach 0.25 ind/km2 24. Coincidentally, this highest density estimate of the N-tiger cat is the same as the lowest ocelot density estimate for Amazonian forests24,29.Tiger cats and margays show high similarity, making misidentifications relatively common39. However, the evaluation of  > 3000 camera trap images of small-medium felids in the Amazon revealed that only one mildly resembled a tiger cat, a finding that supports the species being absent there and does not represent a case of mistaken identity with margays or even ocelots7.The Amazonian range of L. tigrinus is very limited, and populations are expected to be very small. With the upcoming split of L. t. tigrinus and L. t. pardinoides into two different species40, this situation would have serious implications for the conservation of the former. Thus, L. t. tigrinus conservation lies outside the “Amazonian safe haven” of most other carnivore species found there7. The Brazilian drylands Cerrado and Caatinga represent such places for L. t. tigrinus populations. Unfortunately, these biomes have had  > 50% of their cover completely removed41. Very importantly, besides being extremely rare in the Amazonian savannas, this rather limited vegetative formation is also considered highly threatened and of conservation priority42. Therefore, the tiger cat could become an emblematic flagship species representing the uniqueness of this vegetative formation in dire need of protection.In short, the picture that emerges is that although the N-tiger cat uses both rainforests and deciduous forests in the Amazon, it seems to be mostly associated with savanna formations and that its distribution in the Amazon is highly influenced by the ocelot, the dominant mesopredator. The N-tiger cat’s inherent rarity, expected population size, and restricted range in the Amazon suggest that this biome does not in fact represent a safe haven for the global conservation of this small felid. In addition to shedding light on and refining the N-tiger cat distribution in the Amazon, this paper highlights the importance of including biological variables, such as the potential impacts of competitors and predators on species presence and distribution, in SDMs. More

  • in

    A database of seed plants on taxonomy, geography and ecology in the Qinling-Daba Mountains and adjacent areas

    Each of the 23 key variables can be used for analysis. To validate the dataset, we used five plant-related variables (diversity of order, family, genus, species and species endemic to China) to demonstrate the process of using the dataset for analysis as follows:(1) For the four variables of plant taxa “order”, “family”, “genus” and “species”, the similarity and difference in spatial distribution pattern of diversity of different taxa in the Qinling-Daba Mountains climate transition zone were analyzed. The spatial distribution pattern of the diversity of the four taxa is shown in Fig. 3, which is increasingly lower from south (low latitude) to north (high latitude). This result is consistent with the classical latitudinal gradient model of plant diversity. The boundary between higher diversity in the south and lower diversity in the north is roughly located in the area of Funiu Mountains in the eastern Qinling-Daba Mountains, Taibai Mountains in the central Qinling-Daba Mountains and Baishui River in the western Qinling-Daba Mountains. However, with the reduction in taxon scale, the spatial distribution pattern of diversity tends to be complex. Orders (Fig. 3a) and families (Fig. 3b) can be divided by lines, while genera (Fig. 3c) need thicker lines, and species (Fig. 3d) can only be divided by polygons. Figure 3 shows that the taxonomic groups of families are more clearly divided, while species can only be divided by staggered bands. Therefore, when dividing the north–south boundary, the family taxon scale is appropriate, whereas the species scale is more appropriate when studying the north–south transition zone.Fig. 3Spatial distribution of diversity of orders, families, genera and species. (a) The blue dotted line is basically the dividing line of the order diversity of 50 species. The order diversity to the north of the blue dotted line is lower than 50 species, and the order diversity to the south of the blue dotted line is higher than 50 species. (b) The blue dotted line is basically the dividing line of the family diversity of 150 species. The family diversity to the north of the blue dotted line is lower than 150 species, and the family diversity to the south of the blue dotted line is higher than 150 species. (c) The thicker blue dotted line is basically the dividing line of genus diversity of 578–681 species. The genus diversity to the north of the blue dotted line is lower than 578 species, and the genus diversity to the south of the blue dotted line is higher than 681 species. (d) The blue area is basically the dividing line of species diversity of 1385–1618 species. The species diversity to the north of the blue dotted line is lower than 1385 species, and the species diversity to the south of the blue dotted line is higher than 1618 species.Full size imageThe dataset can also count the orders, families and genera that appear in 58 nature reserves, indicating that these orders, families and genera are widely distributed in this area, while the orders, families and genera that only appear in a single nature reserve indicate that these taxa are unique to this nature reserve in this area, reflecting their locality and uniqueness, which is helpful to understanding the specific distribution of plants in detail. The relevant statistics are as follows:
    There are 28 orders present in every nature reserve:
    Liliales, Dipsacales, Lamiales, Fabales, Ericales, Poales, Saxifragales, Malpighiales, Malvales, Asterales, Fagales, Gentianales, Geraniales, Ranunculales, Rosales, Solanales, Apiales, Cornales, Brassicales, Caryophyllales, Dioscoreales, Santalales, Myrtales, Asparagales, Celastrales, Sapindales, Alismatales, and Boraginales.The order that only appears in one nature reserve is Petrosaviales, which appears in the Dabashan Nature Reserve in Chongqing.
    There are 51 families present in every nature reserve:
    Liliaceae, Primulaceae, Plantaginaceae, Lamiaceae, Euphorbiaceae, Cannabaceae, Juncaceae, Fabaceae, Poaceae, Elaeagnaceae, Betulaceae, Apocynaceae, Violaceae, Malvaceae, Crassulaceae, Campanulaceae, Asteraceae, Orchidaceae, Polygonaceae, Orobanchaceae, Onagraceae, Gentianaceae, Geraniaceae, Ranunculaceae, Rubiaceae, Rosaceae, Caprifoliaceae, Thymelaeaceae, Apiaceae, Cyperaceae, Cornaceae, Paeoniaceae, Brassicaceae, Amaryllidaceae, Caryophyllaceae, Rhamnaceae, Santalaceae, Asparagaceae, Celastraceae, Sapindaceae, Adoxaceae, Araliaceae, Berberidaceae, Hydrangeaceae, Scrophulariaceae, Convolvulaceae, Urticaceae, Salicaceae, Papaveraceae, Iridaceae, and Boraginaceae.There are 15 families that only appear in one nature reserve, as shown in Table 2.Table 2 Endemic families of the nature reserves in the Qinling-Daba Mountains and surrounding areas.Full size table
    There are 54 genera present in every nature reserve:
    Patrinia, Polygonum, Sanicula, Plantago, Allium, Delphinium, Euphorbia, Juncus, Cynanchum, Trigonotis, Artemisia, Sorbus, Polygonatum, Scutellaria, Cirsium, Viburnum, Ajuga, Viola, Galium, Geranium, Salix, Epilobium, Gentiana, Ranunculus, Malus, Acer, Rubia, Rosa, Torilis, Lonicera, Adenophora, Philadelphus, Cornus, Paeonia, Rhamnus, Rumex, Carex, Thalictrum, Asparagus, Carpesium, Clematis, Potentilla, Euonymus, Eleutherococcus, Berberis, Spiraea, Rubus, Populus, Vicia, Silene, Iris, Poa, Aster, and Buddleja.There were 225 genera that only appeared in one nature reserve, as shown in Figshare file 269.(2) For the “species endemic to China” variable of plants, we can see from the diversity distribution pattern of species endemic to China in this region (Fig. 4) that the number of endemic species in the Qinling-Daba Mountains is higher than that of species outside of the region, which reflects the strong transition zone in the Qinling-Daba Mountains. The variables of species endemic to China obtained from the Qinling-Daba Mountains and their surroundings were clustered by the Bray–Curtis dissimilarity measure70 and Ward’s minimum variance (the clustering method recommended for plant cluster analysis). The clustering results are shown in Fig. 5a. At the same time, the clustering results are displayed in space. Figure 5b shows that category 3 extends from the east outside the Qinling-Daba Mountains to the Baishuijiang Nature Reserve inside the western Qinling-Daba Mountains, which is consistent with the fact that the Qinling-Daba Mountains are an important ecogeographical “corridor” connecting the east and the west.Fig. 4Spatial distribution of diversity of species endemic to China in the Qinling-Daba Mountains and adjacent areas.Full size imageFig. 5(a) Clustering results of Ward’s connection aggregation of species endemic to China in 58 nature reserves. (b) Spatial distribution of clustering results of species endemic to China; the larger the dot and the darker the color, the earlier it is merged into this category, and the smaller the dot and the lighter the color, the later it is merged into this category.Full size image More

  • in

    Mountain- and brown hare genetic polymorphisms to survey local adaptations and conservation status of the heath hare (Lepus timidus sylvaticus, Nilsson 1831)

    Angerbjörn, A. & Flux, J. E. C. Lepus timidus. Mamm. Species 1–11, https://doi.org/10.2307/3504302 (1995).Bergengren, A. On genetics, evolution and history of distribution of the heath-hare, a distinct population of the Arctic hare, Lepus timidus Lin. Swed. Wildl. (Viltrevy) 6, 381–460 (1969).
    Google Scholar 
    Thulin, C.-G. The distribution of mountain hares Lepus timidus in Europe: a challenge from brown hares L. europaeus? Mamm. Rev. 33, 29–42 (2003).Article 

    Google Scholar 
    Mills, L. S. et al. Camouflage mismatch in seasonal coat color due to decreased snow duration. Proc. Nat.Acad. Sci. 110, 7360–7365 (2013).Article 
    ADS 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Zimova, M. et al. Lack of phenological shift leads to increased camouflage mismatch in mountain hares. Proc.Royal Soc. B: Biol. Sci. 287, 20201786 (2020).Article 

    Google Scholar 
    Levänen, R., Kunnasranta, M. & Pohjoismäki, J. Mitochondrial DNA introgression at the northern edge of the brown hare (Lepus europaeus) range. Ann Zool Fennici 55, 15–24 (2018).Article 

    Google Scholar 
    Thulin, C.-G., Winiger, A., Tallian, A. G. & Kindberg, J. Hunting harvest data in Sweden indicate precipitous decline in the native mountain hare subspecies Lepus timidus sylvaticus (heath hare). J. Nat. Conserv. 64, 126069 (2021).Article 

    Google Scholar 
    Thulin, C.-G., Jaarola, M. & Tegelström, H. The occurrence of mountain hare mitochondrial DNA in wild brown hares. Mol. Ecol. 6, 463–467 (1997).Article 
    CAS 
    PubMed 

    Google Scholar 
    Pohjoismäki, J. L. O., Michell, C., Levänen, R. & Smith, S. Hybridization with mountain hares increases the functional allelic repertoire in brown hares. Sci. Rep. 11, 15771 (2021).Article 
    ADS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Hoekstra, H. E. Genetics, development and evolution of adaptive pigmentation in vertebrates. Heredity (Edinb) 97, 222–234 (2006).Article 
    CAS 

    Google Scholar 
    Hamill, R. M., Doyle, D. & Duke, E. J. Spatial patterns of genetic diversity across European subspecies of the mountain hare, Lepus timidus L. Heredity (Edinb) 97, 355–365 (2006).Article 
    CAS 

    Google Scholar 
    Leach, K., Montgomery, W. I. & Reid, N. Biogeography, macroecology and species’ traits mediate competitive interactions in the order Lagomorpha. Mamm. Rev. 45, 88–102 (2015).Article 

    Google Scholar 
    Marques, J. P. et al. Data Descriptor: Mountain hare transcriptome and diagnostic markers as resources to monitor hybridization with European hares. Sci. Data 4, 1–11 (2017).Article 

    Google Scholar 
    NCBI Sequence Read Archive https://identifiers.org/insdc.sra:SRP358660 (2022).Andrews, S. FastQC: a quality control tool for high throughput sequence data. Babraham Bioinformatics. Preprint at http://www.bioinformatics.babraham.ac.uk/projects/fastqc (2010).Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17, 10–12 (2011).Article 

    Google Scholar 
    Marques, J. P. et al. An annotated draft genome of the mountain hare (Lepus timidus). Genome Biol. Evol. 12, 3656–3662 (2020).Article 
    CAS 
    PubMed 

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

    Google Scholar 
    Broad Institute. Picard toolkit. Broad Institute, GitHub repository. Preprint at https://broadinstitute.github.io/picard/ (2019).Garrison, E. & Marth, G. Haplotype-based variant detection from short-read sequencing. arXiv 1207.3907 (2012).Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156–2158 (2011).Article 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Browning, S. R. & Browning, B. L. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. Am. J. Hum. Genet. 81, 1084–1097 (2007).Article 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Michell, C. T., Pohjoismäki, J. L. O., Spong, G. & Thulin, C.-G. Mountain- and brown hare genetic polymorphisms to survey local adaptations and conservation status of the heath hare (Lepus timidus sylvaticus, Nilsson 1831), Dryad, https://doi.org/10.5061/dryad.3bk3j9kmp (2022).Khan, A. & Mathelier, A. Intervene: a tool for intersection and visualization of multiple gene or genomic region sets. BMC Bioinformatics 18, 287 (2017).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org/.Jombart, T. & Ahmed, I. adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. Bioinformatics 27, 3070–3071 (2011).Article 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Jombart, T. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24, 1403–1405 (2008).Article 
    CAS 
    PubMed 

    Google Scholar 
    Dierckxsens, N., Mardulyn, P. & Smits, G. NOVOPlasty: De novo assembly of organelle genomes from whole genome data. Nucleic Acids Res. 45 (2017).Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Mol. Biol. Evol. 30 (2013).Trifinopoulos, J., Nguyen, L. T., von Haeseler, A. & Minh, B. Q. W-IQ-TREE: a fast online phylogenetic tool for maximum likelihood analysis. Nucleic Acids Res. 44 (2016).Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A. & Jermiin, L. S. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat. Methods 14, 587–589 (2017).Article 
    CAS 
    PubMed 
    PubMed Central 

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

    Google Scholar 
    Paradis, E. & Schliep, K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528 (2019).Article 
    CAS 
    PubMed 

    Google Scholar 
    Ewels, P., Magnusson, M., Lundin, S. & Käller, M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32, 3047–3048 (2016).Article 
    CAS 
    PubMed 
    PubMed Central 

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

    Google Scholar 
    Li, H. Minimap2: Pairwise alignment for nucleotide sequences. Bioinformatics 34, 3094–3100 (2018).Article 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842 (2010).Article 
    CAS 
    PubMed 
    PubMed Central 

    Google Scholar 
    Levänen, R., Thulin, C.-G., Spong, G. & Pohjoismäki, J. L. O. Widespread introgression of mountain hare genes into Fennoscandian brown hare populations. PloS One 13, e0191790 (2018).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Giska, I. et al. The evolutionary pathways for local adaptation in mountain hares. Mol. Ecol. 31, 1487–1503 (2022).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Thulin, C.-G., Isaksson, M. & Tegelström, H. The origin of Scandinavian mountain hares (Lepus timidus). Gibier Faune Savage/Game and Wildlife 14, 463–475 (1997).
    Google Scholar 
    Ferreira, M. S. et al. The legacy of recurrent introgression during the radiation of hares. Syst. Biol. 70, 593–607 (2021).Article 
    PubMed 

    Google Scholar  More

  • in

    Zebras of all stripes repel biting flies at close range

    The evolutionary origins of zebra stripes have been investigated—and debated—for centuries. The trait is rare, conspicuous, and intensely expressed, and thus appears to beg an adaptationist explanation. However, the utility of a complete coat of densely packed, starkly contrasting black-and-white stripes is not immediately apparent. Unlike many conspicuous visual traits, striped pelage is expressed with comparable intensity in both sexes and is thus unlikely to have arisen through sexual selection alone (although in plains zebras, Equus quagga, males have stripes closer to true black than females). Stripes are clearly not aposematic warning signals, nor do they provide camouflage in either the woodland or savannah habitats common across zebra ranges1,2. So, striping presents an ideal evolutionary puzzle: a trait so refined it seems it must be “for” something, but one that confers no clear advantage upon its bearers and imposes apparent costs (conspicuousness) that cannot be explained in Zahavian terms.Scientists have proposed and investigated several possible explanations for the evolution of zebra stripes (reviewed in3). The hypotheses suggest various ways in which stripes may provide a social function (species or individual recognition or social cohesion1,4), a temperature-regulation benefit5,6, an anti-predator effect7,8, or an anti-parasite effect9,10. There is continued debate over both the merits of individual hypotheses and the likelihood of stripes having arisen via a single driver vs. a confluence or alternation of multiple selective pressures6,11.The present study addresses the hypothesis that has thus far received the most empirical support: the anti-parasite hypothesis (also known as the ectoparasite hypothesis12). Zebras, like most ungulates, are harassed by tabanid, glossinid and Stomoxys species of biting flies, which can inflict significant blood loss, transmit disease, and weaken hosts when fly-avoidance behaviors reduce the host’s feeding rate9,13,14. Yet zebras are attacked far less than sympatric ungulates across their African range15,16, and also less than other equids9,17. Zebras also produce odors that may augment their anti-fly defenses18, but so do other sympatric ungulate species18,19, and a host of observations and experiments have demonstrated that black-and-white stripes alone are unattractive, or actively repellent to tabanid, glossinid, and Stomoxys flies17,20,21,22,23.Though the effect of stripes on flies is well-established, the source of the effect remains unexplained. Since Waage’s foundational studies in the 1970s and 1980s9,24 most hypotheses have suggested ways that stripes might interfere with the visual and navigational systems of flies, making it harder for them to locate, identify, or successfully land on striped targets. These hypothetical mechanisms can be roughly grouped by the distance (and the attendant phase of a fly’s orientation and landing behavior) at which they would likely operate:

    From afar: stripes might make it harder for flies to locate and distinguish zebras from background vegetation, perhaps by breaking up their outline9 or varying the way they polarize or reflect light17,31 especially from distances at which composite eyes support only low-resolution vision and cannot resolve zebra stripes as clear bands of alternating color on a single host (estimated at  > 2.0 m22,  > 4.4 m24, and even  > 20 m25).

    At close range (estimates range from 0.5 to 4.0 m26): stripes might interfere with orientation or landing behavior via any of several disruptive or ‘dazzle’-related visual effects27. For example, stripes might affect ‘optic flow’, or the fly’s perceived relative motion to its target as it approaches, by creating an illusion of false direction or speed of motion (e.g., via variants of the ‘barber pole’ or ‘wagon wheel’ effects28). Alternatively, relative motion to a striped pattern within the visual field may create the perception of self-rotation, inducing the fly’s involuntary ‘optomotor response’ and resulting in an avoidance turn in an effort to stay on a straight course29.

    Finally, stripes might cause confusion in the transition between long- and short-distance orientation. If zebras appear as blurred gray from a distance and then, at closer range, suddenly resolve into a sequence of floating black and white bars, this abrupt ‘visual transformation’26 might disrupt the behavioral sequence that facilitates landing.

    Within these categories, hypotheses have proliferated faster than experimental tests of many of the proposed mechanisms. The very active literature on this question has grown in somewhat haphazard fashion, as curious researchers test new possibilities without eliminating old ones6. Importantly, few experiments have controlled the distance from which flies are first able to view potential landing sites (but see23). While growing evidence supports a mechanism operating at close range22,26, failing to restrict the starting distance of the fly means that the full set of possible mechanisms outlined above all remain plausible contributors to most previous results.Additionally, while many studies have, appropriately, used artificial stimuli to isolate basic effects of color, pattern, brightness, and light polarization of (usually flat) test surfaces, possible contributions of several aspects of natural zebra pelage remain untested. Controlled experiments have used various landing substrates, including striped and solid oil tray traps, sticky plastic, smooth plastic17, cloth (Experiment 2 in22), horse blankets or sheets26, and paint on live animals30. These have all clearly demonstrated a broadly replicable visual effect: stripes, and some other juxtapositions of black and white (e.g., checkerboard patterns26), repel flies. However, insofar as specific features of zebra pelage factor into proposed mechanisms of fly repellence—the reflective properties of “smooth, shiny” coats31; the orientation of the stripes17,32; the light-polarizing effects of black and white hair vs. background vegetation25; and the complex structure of hair25—there is a need for more experiments that present natural targets to wild flies (but see22,33). Similarly, most experiments have compared landing preferences between black-and-white striped, solid black, solid white, and occasionally solid grey substrates, which have served as important controls for determining that light polarization, rather than a combination of polarization and brightness, is sufficient to induce the effect of stripe avoidance17. However, it is now time to refocus on the original question by presenting flies with more realistic choices. Since biting flies seeking a bloodmeal on the African savannah seldom encounter solid black hosts, and even more rarely solid white hosts, landing choices should be compared between zebra stripes and common coat colors of sympatric mammals, namely various shades of brown. Further, tabanid, glossinid, and Stomoxys flies all avoid landing on stripes that are the same width or narrower than the widest zebra stripes 17,23, and there is some evidence that narrower stripes are even more repellent to tabanids17. This pattern is potentially significant in the application of the anti-parasite hypothesis to an adaptive explanation for the striking variation in stripe width across zebra species and between the different areas of the body on individual zebras22, but must first be confirmed with experiments using real zebra pelage.Here, we present a simple experiment designed to address each of these gaps in the literature on the anti-fly benefits of zebra stripes. In this field experiment, the landing choices of flies were tested entirely within the range at which all estimates agree flies should be able to perceive the presented stripes ( More

  • in

    Low levels of sibship encourage use of larvae in western Atlantic bluefin tuna abundance estimation by close-kin mark-recapture

    Our results show that GoM BFT larval survey samples can provide the crucial mark events for eventual CKMR estimates of adult abundance. The adult parents marked by larval samples can be directly recaptured in the fishery many years later as POPs, and indirectly through their progeny in future samples of larvae, as evidenced by the two cross-cohort HSPs (XHSPs) recovered in this study, which imply that a parent survived and spawned in the GoM in consecutive years. As more cohorts are sampled in future, the growing number of XHSPs could be used to estimate average adult survival rates, in addition to helping with the estimation of adult abundance31, as is now done for southern blue tuna40.There is a modest level of sibship within our 2016 samples, and a high level (involving over half the samples) in 2017, but it turns out not to be high enough to cause serious problems for POP-based CKMR. High sibship per se does not lead to bias in CKMR by virtue of the statistical construction of the estimate, but it does increase variance, which can be summarized through a reduction in effective sample size. In a POP-based CKMR model, our effective sample size would be about 75% of nominal for the two years combined, or 66% of nominal for the targeted sampling of 2017. Since it is actually the product of adult and juvenile sample sizes which drives precision in CKMR14, one way to think about the 75% is that we will need about 33% more adult samples to achieve a given precision on abundance estimates than if we had somehow been able to collect the same number of “independent” juvenile samples (i.e. without oversampling siblings). That increase is appreciable but entirely achievable; for WBFT, it is logistically much easier to collect more feeding-ground adult samples than to collect more larvae, and at present there is no known practical way to collect large numbers of older, more dispersed, and thus more independent, juvenile western origin bluefin tuna (WBFT).This study was motivated by the concern that sibship might be a serious impediment to use of WBFT larvae for CKMR. High levels of sibship have been found in larval collections for other taxa despite a pelagic larval phase, suggesting that abiotic factors can impede random mixing of larvae after a spawning event41. Our larval samples were only a few days old (4–11) and thus had little time to disperse since fertilization; our concern beforehand was that each tow might sample the offspring of a very small number of adults (one spawning group in one night), and in 2017 that repeatedly towing the same water mass might simply be resampling the same “family”. In practice, though, the cumulative effect was limited. Samples were not dominated by progeny from just a few adults; the maximum DPG size (i.e., number of offspring from any one adult) was 5, which is under 2% of the larval sample size. There are several possible reasons for this finding. First, plankton sample tows are typically standardized to a ten-minute duration, covering on average about 0.3 nautical miles. Based on continuous plankton cameras42, each tow is likely to tow through multiple patches of zooplankton, and therefore potentially multiple patches of BFT larvae. Second, spawning aggregations of BFT may contain many adults. For example, on the spawning grounds near the Balearic Islands in the Mediterranean, purse seine fisheries target spawning fish and individual net sets routinely capture upwards of 500 mature individuals43. These numbers suggest that BFT spawner aggregations can be quite large, although the number of individuals that contribute gametes to a single spawning event may be lower. The results of this study pose intriguing scenarios for understanding BFT larval ecology and spawning behavior, which could be explored with larger sample sizes paired with data on oceanographic conditions, direct observation of spawning aggregations, and modeling to compare observed and predicted dispersal. The results of this study are based on just two years of sampling, and numerous practical and theoretical challenges remain to fully understand BFT reproduction in the GoM.Our sibship impact calculations assume use of an unmodified adult-size-based CKMR POP model, where each juvenile is compared to each adult taking into account the latter’s size (e.g.,14). That will give unbiased estimates, which we regard as essential in a CKMR model. However, for WBFT the estimates are not fully statistically efficient, in that some adults receive more statistical weight than others because they are marked more often (by having a large DPG), and thus variance might not be the lowest achievable. Modifying the model to fix that would be simple in a “cartoon” CKMR setting where all adults are identical (e.g., Fig. 1 of14), simply by first condensing each DPG to a single representative, then only using those representatives (rather than all the larvae) in POP comparisons. Each marked parent then receives the same weight, giving maximum efficiency. For the cartoon, this condensed-DPG model still gives an unbiased estimate of abundance, because each DPG has one parent of given sex, and the chance of any sampled cartoon adult of that sex being that parent is 1/N. The DPG-condensed effective sample size is simply half the number of distinct parents, which would be a little larger than the effective sample sizes for the unmodified model shown in Table 3; e.g., in 2017, 504/2 = 252 versus 209. However, no such straightforward improvement is available for an adult-size-based CKMR model such as is needed for WABFT. Using condensed DPGs directly would bias the juvenile sampling against larger more-fecund adults, whose DPGs will tend on average to be larger and thus to experience disproportionate condensation. Those adults would be marked less often by the DPG-condensed juveniles than the model assumes, violating the basic requirements for unbiased CKMR in14. A more sophisticated model might be able to combine unbiasedness with higher efficiency but, since the unmodified adult-size-based POP model that we expect to use is unbiased and only mildly inefficient (at worst 209/252 = 83% efficient, in 2017) there seems no particular need for extra complications at present. However, that may not hold true if we eventually move to a POP + XHSP model, where the impact on unmodified CKMR variance is worse (though there is still no bias, for the same reason as with POPs). Intuitively, the biggest impact that a DPG of size 5 can have in a POP model is to suddenly raise the number of POPs by 5 if its parent happens to be sampled; within a useful total of, say, 75 POPs, the influence is not that large. But if two DPGs both of size 5 in different cohorts happen to share a parent, then the total of XHSPs suddenly jumps by 25— likely a substantial proportion of total XHSPs. Supplementary Material B also includes effective sample size formulae for a simplified XHSP-only model, which demonstrate the increased impact of within-cohort sibship; for our WBFT samples, it turns out that the XHSP-effective size is slightly lower for the targeted 2017 samples (110) than for the 2016 samples (130), unlike the POP-only effective size. Dropping from a maximum theoretical effective sample size of 252 (half the number of DPGs) down to 110 would be rather inefficient and would increase the number of years of sampling required to yield a useful XHSP dataset. This motivates developing a modified POP + XHSP model that retains unbiasedness without sacrificing too much efficiency. In principle, that can be done by condensing each DPG but then conditioning its comparison probabilities on the DPG’s original size, in accordance with the framework in14. This is a topic for subsequent research, and the results will inform future sampling strategy decisions for WBFT.One potential difficulty for western BFT CKMR might occur if a substantial proportion of animals reaching maturity are the offspring of “Western” (in genetic terms) adults who persistently spawn in the western North Atlantic but outside the GoM. However, as long as the adults marked by GoM larvae are well mixed at the time of sampling with any western adults that do spawn outside of the GoM, the total POP-based population estimate of genetically-western BFT from CKMR will remain unbiased. Given evidence from tagging of widespread adult movements within the western North Atlantic2, good mixing in the sampled feeding grounds seems likely; so, even if successful non-GoM western BFT spawning really is commonplace, there should not be a problem with relying on GoM larvae for at least the POP component of CKMR14.Studies of fish early life history have long been considered to have great potential to provide novel insight into the unique population dynamics of fishes44,45,46. Sampling efforts aimed at estimating fish recruitment dynamics have spawned a diversity of larval survey programs. Examples of these long-term programs include the California Cooperative Oceanic Fisheries Investigations, International Council for the Exploration of the Sea (ICES) surveys in the North Atlantic and adjacent areas, Southeast Monitoring and Assessment Program (SEAMAP) in the GoM, Ecosystem Monitoring (EcoMon) in the Northeast U.S., and numerous others, many of which provide indices of larval abundance widely used in fisheries and ecosystem assessments. Yet, as a result of the inherent patchiness of larvae42, sampling variability, and highly variable density dependent mortality45, fisheries scientists have often struggled to determine how larval surveys relate to the adult fish populations. Inclusion of estimates of sibship among larvae collected in surveys could refine estimates of adult spawning stock biomass estimated from these surveys.The results of this study also represent products of decades of work and coordination in obtaining high-quality DNA from larval specimens. Key steps to successful genotyping of larvae include ensuring that larvae are preserved, sorted, and handled in 95% non-denatured ethanol. In addition, strict instrument cleaning protocols must be followed, and stomachs should be removed or avoided (this study used larval tails and, when possible, eyes to avoid cross contamination of prey contents, including possible congeners and other BFT individuals). Exposure to hot lamps during the sorting and dissection processes should also be minimized to ensure that DNA quality is sufficiently high for genotyping-by-sequencing. Although the tissues available for genetic analysis were limited by the needs of other experiments that required BFT tissues, otoliths, gut contents, and other information from the same larvae, we were able to successfully genotype most larvae greater than 6 mm SL and identify thousands of informative SNPs. The lower size limit of larvae could likely be decreased if whole specimens were available for genotyping, although the use of younger larvae could increase the incidence of sibship.In summary, while we observed both FSPs and HSPs in larval collections, with elevated sibship overall and with siblings being more prevalent within tows and in nearby tows, the level of sibship was sufficiently low that collections of GoM BFT larvae can still provide the critical genetic mark of parental genotypes required for CKMR. Our results demonstrate a crucial proof of concept and are the first step towards an operational CKMR modelling estimate of spawning stock abundance for western BFT. More

  • in

    Single-cell measurements and modelling reveal substantial organic carbon acquisition by Prochlorococcus

    Isotope labelling and phylogenetic analysis of a natural marine bacterioplankton population at seaMediterranean seawater was collected during August 2017 (station N1200, 32.45° N, 34.37 °E) from 11 depths by Niskin bottles and divided into triplicate 250 ml polycarbonate bottles. Two bottles from each depth were labelled with 1 mM sodium bicarbonate-13C and 1 mM ammonium-15N chloride (Sigma-Aldrich), and all three bottles (two labelled and one control) were incubated at the original depth and station at sea for 3.5 h around mid-day. The stable isotopes were chosen to enable direct comparison of C and N uptake in single cells, and the short incubation time was chosen to minimize isotope dilution and potential recycling and transfer of 13C and 15N between community members25. After incubation, bottles were brought back on board and the incubations were stopped by fixing with 2× electron-microscopy-grade glutaraldehyde (2.5% final concentration) and stored at 4 °C until sorting analysis. Cell sorting, NanoSIMS analyses and the calculation of uptake rates were performed as described in Roth-Rosenberg et al.26.DNA collection and extraction from seawaterSamples for DNA were collected on 0.22 µm Sterivex filters (Millipore). Excess water was removed using a syringe, 1 ml lysis buffer (40 mM EDTA, 50 mM Tris pH 8.3, and 0.75 M sucrose) was added and both ends of the filter were closed with parafilm. Samples were kept at −80 °C until extraction. DNA was extracted by using a semi-automated protocol including manual chemical cell lysis before automated steps using the QIAamp DNA Mini Protocol: DNA Purification from Blood or Body Fluids (Spin Protocol, starting from step 6, at the BioRap unit, Faculty of Medicine, Technion). The manual protocol began with thawing the samples, then the storage buffer was removed using a syringe and 170 µl lysis buffer added to the filters. Thirty microlitres of Lysozyme (20 mg ml−1) were added to the filters and incubated at 37 °C for 30 min. After incubation, 20 µl proteinase K and 200 µl buffer AL (from the Qiagen kit) were added to the tube for 1 h at 56 °C (with agitation). The supernatant was transferred to a new tube, and DNA was extracted using the QIAcube automated system. All DNA samples were eluted in 100 μl DNA-free distilled water.ITS PCR amplificationPCR amplification of the ITS was carried out with specific primers for Prochlorococcus CS1_16S_1247F (5′-ACACTGACGACATGGTTCTACACGTACTACAATGCTACGG) and Cs2_ITS_Ar (5′-TACGGTAGCAGAGACTTGGTCTGGACCTCACCCTTATCAGGG)21,22. The first PCR was performed in triplicate in a total volume of 25 μl containing 0.5 ng of template, 12.5 μl of MyTaq Red Mix (Bioline) and 0.5 μl of 10 μM of each primer. The amplification conditions comprised steps at 95 °C for 5 min, 28/25 (16 S/ITS) cycles at 95 °C for 30 s, 50 °C for 30 s and 72 °C for 1 min followed by one step of 5 min at 72 °C. All PCR products were validated on a 1% agarose gel, and triplicates were pooled. Subsequently, a second PCR amplification was performed to prepare libraries. These were pooled and after a quality control sequenced (2 × 250 paired-end reads) using an Illumina MiSeq sequencer. Library preparation and pooling were performed at the DNA Services facility, Research Resources Center, University of Illinois at Chicago. MiSeq sequencing was performed at the W.M. Keck Center for Comparative and Functional Genomics at the University of Illinois at Urbana-Champaign.ITS sequence processingPaired-end reads were analysed using the Dada2 pipeline46. The quality of the sequences per sample was examined using the Dada2 ‘plotQualityProfile’ command. Quality filtering was performed using the Dada2 ‘filterAndTrim’ command with parameters for quality filtering truncLen=c(290,260), maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE, trimLeft=c(20,20). Following error estimation and dereplication, the Dada2 algorithm was used to correct sequences. Merging of the forward and reverse reads was done with minimum overlap of 4 bp. Detection and removal of suspected chimaeras was done with command ‘removeBimeraDenovo’. In total, 388,417 sequences in 484 amplicon sequence variants were counted. The amplicon sequence variants were aligned in MEGA6 (ref. 47), and the first ~295 nucleotides, corresponding to the 16S gene, were trimmed. The ITS sequences were then classified using BLASTn against a custom database of ITS sequences from cultured Prochlorococcus and Synechococcus strains as well as from uncultured HL and LL clades.Individual-based modelPlanktonIndividuals.jl (v0.1.9) was used to run the individual-based simulations48. Briefly, the cells fix inorganic carbon through photosynthesis and nitrogen, phosphorus and DOC from the water column into intracellular quotas and grow until division or grazing. Cell division is modelled as a probabilistic function of cell size. Grazing is represented by a quadratic probabilistic function of cell population. Cells consume nutrient resources, which are represented as Eulerian, density-based tracers. A full documentation of state variables and model equations are available online at https://juliaocean.github.io/PlanktonIndividuals.jl/dev/. Equations related to mixotrophy are shown below as an addition to the online documentation.$$V_{{mathrm{DOC}}} = V_{{mathrm{DOC}}}^{{mathrm{max}}} cdot {{mathrm{max}}}left( {0.0,{{mathrm{min}}}left( {1.0,,frac{{q_{mathrm{C}}^{{mathrm{max}}} – q_{mathrm{C}}}}{{q_{mathrm{C}}^{{mathrm{max}}} – q_{mathrm{C}}^{{mathrm{min}}}}}} right)} right) cdot frac{{{mathrm{DOC}}}}{{{mathrm{DOC}} + K_{{mathrm{DOC}}}^{{mathrm{sat}}}}}$$
    (1)
    $$f_{{mathrm{PS}}} = frac{{P_{mathrm{S}}}}{{P_{mathrm{S}} + V_{{mathrm{DOC}}}}}$$
    (2)
    $$V_{{mathrm{DOC}}} = 0,,{mathrm{if}},f_{{mathrm{PS}}} < f_{{mathrm{PS}}}^{{mathrm{min}}}$$ (3) where VDOC is the cell-specific DOC uptake rate (mol C cell−1 s−1), (V_{{mathrm{DOC}}}^{{mathrm{max}}}) is the maximum cell-specific DOC uptake rate (mol C cell−1 s−1), (q_{mathrm{C}}^{{mathrm{max}}}) is the maximum cell carbon quota (mol C cell−1), (q_{mathrm{C}}^{{mathrm{min}}}) is the minimum cell carbon quota (mol C cell−1). The maximum and minimum functions here is used to keep qC between (q_{mathrm{C}}^{{mathrm{min}}}) and (q_{mathrm{C}}^{{mathrm{max}}}). (K_{{mathrm{DOC}}}^{{mathrm{sat}}}) is the half-saturation constant for DOC uptake (mol C m−3). fPS is the fraction of fixed C originating from photosynthesis (PS, mol C cell−1 s−1). DOC uptake stops when fPS is smaller than (f_{{mathrm{PS}}}^{{mathrm{min}}})(minimum fraction of fixed C originating form photosynthesis, 1% by default) according to laboratory studies of Prochlorococcus that showed that they cannot survive long exposure to darkness (beyond several days) even when supplied with organic carbon sources13. (1 − fPS) is also shown in Fig. 3 as the contribution of DOC uptake.We set up two separate simulations; each of them has a population of either an obligate photo-autotroph or a mixotroph that also consumes DOC. The initial conditions and parameters (Supplementary Table 3) are the same for the two simulations except the ability of mixotrophy. The simulations were run with a timestep of 1 min for 360 simulated days to achieve a steady state. We run the two simulations for multiple times in order to get the range of the stochastic processes.Evaluation of autotrophic growth ratesWe evaluated the carbon-specific, daily-averaged carbon fixation rate, ℙ as a function of light intensity (I, µE), following Platt et al.33:$${Bbb P} = frac{1}{{Delta t}}{int}_0^{Delta t} {frac{{q_{{mathrm{Chl}}}}}{{q_{mathrm{C}}}}} P_{mathrm{S}}^{{mathrm{Chl}}}left( {1 - e^{ - alpha _{{mathrm{Chl}}}I/P_{mathrm{S}}^{{mathrm{Chl}}}}} right)e^{ - beta _{{mathrm{Chl}}}I/P_{mathrm{S}}^{{mathrm{Chl}}}}Delta t$$ (4) Here, (P_{mathrm{S}}^{{mathrm{Chl}}}), αChl and βChl are empirically determined coefficients representing the chlorophyll-a-specific carbon fixation rate (mol C (mol Chl)−1 s−1), the initial slope of the photosynthesis–light relationship and photo-inhibition effects at high photon fluxes, respectively. We impose empirically determined values for (P_{mathrm{S}}^{{mathrm{Chl}}}), αChl and βChl from the published study of Moore and Chisholm24. The natural Prochlorococcus community comprises HL and LL ecotypes, which have different values of (P_{mathrm{S}}^{{mathrm{Chl}}}), αChl and βChl, and the community growth rate is expected to be between that of HL extremes and LL extremes. Therefore, we use photo-physiological parameters for an HL-adapted ecotype (MIT9215), acclimated at 70 µmol photons m−2 s−1 and an LL-adapted ecotype (MIT9211), acclimated 9 µmol photons m−2 s−1. The models with these values are shown as the different lines in Fig. 2b,d. I is the hourly PAR, estimated by scaling the observed noon value at each depth with a diurnal variation evaluated from astronomical formulae based on geographic location and time of year37,38.(frac{{q_{{mathrm{Chl}}}}}{{q_{mathrm{C}}}}) is the molar chlorophyll-a to carbon ratio, which is modelled as a function of growth rate and light intensity using the Inomura34 model (equation 17 therein) where parameters were calibrated with laboratory data from Healey49. In addition, the maximum growth rate ((mu _{{mathrm{max}}}^I)) based on macromolecular allocation is also estimated using the Inomura model (equation 30 therein). An initial guess of the growth rate and the empirically informed light intensity are used to estimate (frac{{q_{{mathrm{Chl}}}}}{{q_{mathrm{C}}}}), which is then used to evaluate the light-limited, photoautotrophic growth rate$${Bbb V}_{mathrm{C}}^{{mathrm{auto}}} = min left( {{Bbb P} - K_{mathrm{R}},mu _{{mathrm{max}}}^I} right)$$ (5) from which the (frac{{q_{{mathrm{Chl}}}}}{{q_{mathrm{C}}}}) is again updated. The light-limited growth rate is used to re-evaluate the (frac{{q_{{mathrm{Chl}}}}}{{q_{mathrm{C}}}}). Repeating this sequence until the values converge, ({Bbb V}_{mathrm{C}}^{{mathrm{auto}}}) and (frac{{q_{{mathrm{Chl}}}}}{{q_{mathrm{C}}}}) are solved iteratively.The nitrogen-specific uptake rate of fixed nitrogen (day−1) is modelled as$${Bbb V}_{{{mathrm{N}}}} = {Bbb V}_{mathrm{N}}^{{mathrm{max}}}frac{1}{{Q_{mathrm{N}}}}frac{N}{{N + K_{{{mathrm{N}}}}}}$$ (6) where values of the maximum uptake rate, ({Bbb V}_{mathrm{N}}^{{mathrm{max}}}), and half-saturation, KN, are determined from empirical allometric scalings35, along with a nitrogen cell quota QN from Bertilsson et al.39.The P-limited growth rate, or the phosphorus-specific uptake rate of phosphate (day−1), is modelled as$${Bbb V}_{mathrm{P}} = {Bbb V}_{mathrm{P}}^{{mathrm{max}}}frac{1}{{Q_{mathrm{P}}}}frac{{{mathrm{PO}_{4}}^{3 - }}}{{{mathrm{PO}_{4}}^{3 - } + K_{mathrm{P}}}}$$ (7) where values of the maximum uptake rate, ({Bbb V}_{mathrm{P}}^{{mathrm{max}}}). and half-saturation, KP, are determined from empirical allometric scalings35, along with a nitrogen cell quota QP from Bertilsson et al.39.Iron uptake is modelled as a linear function of cell surface area (SA), with rate constant ((k_{{mathrm{Fe}}}^{{mathrm{SA}}})) following Lis et al.36.$${Bbb V}_{{mathrm{Fe}}} = k_{{mathrm{Fe}}}^{{mathrm{SA}}} cdot {mathrm{SA}}frac{1}{{Q_{{mathrm{Fe}}}}}{mathrm{Fe}}$$ (8) The potential light-, nitrogen-, phosphorus- and iron-limited growth rates (({Bbb V}_{mathrm{C}},{Bbb V}_{mathrm{N}},{Bbb V}_{mathrm{P}},{Bbb V}_{{mathrm{Fe}}})) were evaluated at each depth in the water column and the minimum is the local modelled photo-autotrophic growth rate estimate, assuming no mixotrophy (Fig. 2b,d, blue lines). Parameters used in this evaluation are listed in Supplementary Table 2.An important premise of this study is that heterotrophy is providing for the shortfall in carbon under very low light conditions, but not nitrogen. It is known that Prochlorococcus can assimilate amino acids9 and therefore the stoichiometry of the heterotrophic contribution might alter the interpretations. However, it is also known that Prochlorococcus can exude amino acids40, which might cancel out the effects on the stoichiometry of Prochlorococcus.For the estimates of phototrophic growth rate from local environmental conditions (Fig. 2) we employed photo-physiological parameters from laboratory cultures of Prochlorococcus24. For the purposes of this study, we have assumed that the photosynthetic rates predicted are net primary production, which means that autotrophic respiration has been accounted for in the measurement. However, the incubations in that study were of relatively short timescale (45 min), which might suggest they are perhaps more representative of gross primary production. If this is the case, our estimates of photo-autotrophic would be even lower after accounting for autotrophic respiration, and thus would demand a higher contribution from heterotrophic carbon uptake. In this regard, our estimates might be considered a lower bound for organic carbon assimilation.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More