Philippa Kaur

Ryther, J. H. Photosynthesis and fish production in the sea. Sci. (80-.) 166, 72–76 (1969).ADS 
CAS 
Article 

Google Scholar 
Follows, M. J., Dutkiewicz, S., Grant, S. & Chisholm, S. W. Emergent biogeography of microbial communities in a model ocean. Sci. (80-.). 315, 1843–1846 (2007).ADS 
CAS 
Article 

Google Scholar 
Edwards, K. F., Litchman, E. & Klausmeier, C. A. Functional traits explain phytoplankton community structure and seasonal dynamics in a marine ecosystem. Ecol. Lett. 16, 56–63 (2013).PubMed 
Article 

Google Scholar 
Nemergut, D. R. et al. Patterns and processes of microbial community assembly. Microbiol. Mol. Biol. Rev. 77, 342–356 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
Villarino, E. et al. Large-scale ocean connectivity and planktonic body size. Nat. Commun. 9, 142 (2018).Collins, S., Rost, B. & Rynearson, T. A. Evolutionary potential of marine phytoplankton under ocean acidification. Evol. Appl. 7, 140–155 (2014).CAS 
PubMed 
Article 

Google Scholar 
Rusch, D. B. et al. The Sorcerer II global ocean sampling expedition: Northwest Atlantic through Eastern Tropical Pacific. PLOS Biol. 5, e77 (2007).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
de Vargas, C. et al. Eukaryotic plankton diversity in the sunlit ocean. Sci. (80-.). 348, 1261605–1/11 (2015).Sunagawa, S. et al. Structure and function of the global ocean microbiome. Sci. (80-.) 348, 1–10 (2015).Article 
CAS 

Google Scholar 
Fuhrman, J. A. et al. A latitudinal diversity gradient in planktonic marine bacteria. Proc. Natl Acad. Sci. 105, 7774–7778 (2008).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Righetti, D., Vogt, M., Gruber, N., Psomas, A. & Zimmermann, N. E. Global pattern of phytoplankton diversity driven by temperature and environmental variability. Sci. Adv. 5, 1–11 (2019).Article 

Google Scholar 
Cermeño, P. et al. The role of nutricline depth in regulating the ocean carbon cycle. PNAS 105, 20344–20349 (2008).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Barton, A. D., Dutkiewicz, S., Flierl, G., Bragg, J. & Follows, M. J. Patterns of diversity in marine phytoplankton. Sci. (80-.) 327, 1509–1511 (2010).ADS 
CAS 
Article 

Google Scholar 
Mantyla, A. W., Venrick, E. L. & Hayward, T. L. Primary production and chlorophyll relationships, derived from ten year of CalCOFI measurements. Calif. Cooperative Ocean. Fish. Investig. Rep. 36, 159–166 (1995).
Google Scholar 
Hayward, T. L. & Venrick, E. L. Nearsurface pattern in the California Current: Coupling between physical and biological structure. Deep. Res. Part II Top. Stud. Oceanogr. https://doi.org/10.1016/S0967-0645(98)80010-6 (1998).Article 

Google Scholar 
Venrick, E. L. Floral patterns in the California Current: The coastal-offshore boundary zone. J. Mar. Res. 67, 89–111 (2009).Article 

Google Scholar 
Powell, J. R. & Ohman, M. D. Covariability of zooplankton gradients with glider-detected density fronts in the Southern California Current System. Deep Sea Res. Part II Top. Stud. Oceanogr. 112, 79–90 (2015).ADS 
CAS 
Article 

Google Scholar 
Taylor, A. G., Landry, M. R., Selph, K. E. & Wokuluk, J. J. Temporal and spatial patterns of microbial community biomass and composition in the Southern California Current Ecosystem. Deep. Res. Part II Top. Stud. Oceanogr. 112, 117–128 (2015).Catlett, D. et al. Diagnosing seasonal to multi-decadal phytoplankton group dynamics in a highly productive coastal ecosystem. Prog. Oceanogr. 197, 102637 (2021).Article 

Google Scholar 
Lilly, L. E. & Ohman, M. D. CCE IV: El Niño-related zooplankton variability in the southern California Current System. Deep. Res. Part I Oceanogr. Res. Pap. 140, 36–51 (2018).ADS 
Article 

Google Scholar 
Richardson, A. J. et al. Using continuous plankton recorder data. Prog. Oceanogr. 68, 27–74 (2006).ADS 
Article 

Google Scholar 
Wang, Z. et al. Microbial communities across nearshore to offshore coastal transects are primarily shaped by distance and temperature. Environ. Microbiol. 1462–2920.14734. https://doi.org/10.1111/1462-2920.14734 (2019).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, S161–S179 (2020).Ibarbalz, F. M. et al. Global Trends in Marine Plankton Diversity across Kingdoms of Life. Cell 1084–1097. https://doi.org/10.1016/j.cell.2019.10.008 (2019).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 
Gilbert, J. A. et al. Defining seasonal marine microbial community dynamics. ISME J. 6, 298–308 (2012).CAS 
PubMed 
Article 

Google Scholar 
Karl, D. M. & Lukas, R. The Hawaii Ocean Time-series (HOT) program: background, rationale and field implementation. Deep. Res. Part II Top. Stud. Oceanogr. 43, 129–156 (1996).ADS 
CAS 
Article 

Google Scholar 
Steinberg, D. K. et al. Overview of the US JGOFS Bermuda Atlantic Time-series Study (BATS): A decade-scale look at ocean biology and biogeochemistry Overview of the US JGOFS Bermuda Atlantic Time-series Study (BATS): a decade-scale look at ocean biology and biogeochemistry. Deep. Res. Part II Top. Stud. Oceanogr. 48, 1405–1447 (2015).ADS 
Article 

Google Scholar 
Needham, D. M. & Fuhrman, J. A. Pronounced daily succession of phytoplankton, archaea and bacteria following a spring bloom. Nat. Microbiol. 1, 16005 (2016).Zhu, Z. et al. Understanding the blob bloom: Warming increases toxicity and abundance of the harmful bloom diatom Pseudo-nitzschia in California coastal waters. Harmful Algae 67, 36–43 (2017).CAS 
PubMed 
Article 

Google Scholar 
Mcclatchie, S. et al. State of the California Current 2015–16: Comparisons with the 1997–98 El Niño. Calif. Cooperative Ocean. Fish. Investig. Rep. 57, (2016).Walker, H. J. Jr et al. Unusual occurrences of fishes in the Southern California Current System during the warm water period of 2014–2018. Estuar. Coast. Shelf Sci. 236, 106634 (2020).Article 

Google Scholar 
Kahru, M., Jacox, M. G. & Ohman, M. D. CCE1: Decrease in the frequency of oceanic fronts and surface chlorophyll concentration in the California Current System during the 2014–2016 northeast Pacific warm anomalies. Deep. Res. Part I Oceanogr. Res. Pap. 140, 4–13 (2018).ADS 
Article 

Google Scholar 
Azam, F. et al. The Ecological Role of Water-Column Microbes in the Sea. Mar. Ecol. Prog. Ser. 10, 257–263 (1983).ADS 
Article 

Google Scholar 
Calbet, A. & Landry, M. R. Phytoplankton growth, microzooplankton grazing, and carbon cycling in marine systems. Limnol. Oceanogr. 49, 51–57 (2004).ADS 
CAS 
Article 

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 
Kohonen, T. Exploration of very large databases by self-organizing maps. IEEE Int. Conf. Neural Networks – Conf. Proc. 1, (1997).Istvánovics, V. Eutrophication of Lakes and Reservoirs. Encycl. Inl. Waters 157–165 https://doi.org/10.1016/B978-012370626-3.00141-1 (2009).Partensky, F., Blanchot, J. & Vaulot, D. Differential distribution and ecology of Prochlorococcus and Synechococcus in oceanic waters: a review. Bull. Oceanogr. Monaco 19, 457–475 (1999).
Google Scholar 
Laws, E. A., Falkowski, P. G., Smith, W. O., Ducklow, H. & McCarthy, J. J. Temperature effects on export production in the open ocean. Global Biogeochem. Cycles 14, (2000).Grover, J. P. Resource Competition in a Variable Environment: Phytoplankton Growing According to Monod’s Model. Am. Nat. 136, 771–789 (1990).Article 

Google Scholar 
Benincá, E. et al. Chaos in a long-term experiment with a plankton community. Nature 451, 822–825 (2008).ADS 
PubMed 
Article 
CAS 

Google Scholar 
Williams, R. G. & Follows, M. J. Ocean Dynamics and the Carbon Cycle: Principles and Mechanisms. Book (2011).Lindegren, M., Checkley, D. M., Ohman, M. D., Koslow, J. A. & Goericke, R. Resilience and stability of a pelagic marine ecosystem. Proc. R. Soc. B Biol. Sci. 283, (2016).Vallina, S. M. et al. Global relationship between phytoplankton diversity and productivity in the ocean. Nat. Commun. 1–10 https://doi.org/10.1038/ncomms5299 (2014).Chase, J. M. & Leibold, M. A. Spatial scale dictates the productivity-biodiversity relationship. Nature 416, 427–430 (2002).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Jacox, M. G., Edwards, C. A., Hazen, E. L. & Bograd, S. J. Coastal Upwelling Revisited: Ekman, Bakun, and Improved Upwelling Indices for the U.S. West Coast. J. Geophys. Res. Ocean. 123, 7332–7350 (2018).ADS 
Article 

Google Scholar 
Zaba, K. D. & Rudnick, D. L. The 2014-2015 warming anomaly in the Southern California Current System observed by underwater gliders. Geophys. Res. Lett. 43, 1241–1248 (2016).ADS 
Article 

Google Scholar 
Weber, E. D. et al. State of the California Current 2019–2020: Back to the Future With Marine Heatwaves? Front. Mar. Sci. 8, (2021).Closset, I. et al. Diatom response to alterations in upwelling and nutrient dynamics associated with climate forcing in the California Current System. Limnol. Oceanogr. 1–16. https://doi.org/10.1002/lno.11705 (2021).Kenitz, K. M. et al. Environmental drivers of population variability in colony-forming marine diatoms. Limnol. Oceanogr. 65, 2515–2528 (2020).ADS 
Article 

Google Scholar 
Mullin, M. M. Biomasses of large-celled phytoplankton and their relation to the nitricline and grazing in the California current system off Southern California, 1994–1996. Calif. Cooperative Ocean. Fish. Investig. Rep. 39, 117–123 (1998).
Google Scholar 
Rykaczewski, R. R. & Checkley, D. M. Influence of ocean winds on the pelagic ecosystem in upwelling regions. PNAS 105, 1965–1970 (2007).ADS 
Article 

Google Scholar 
Grzymski, J. J. & Dussaq, A. M. The significance of nitrogen cost minimization in proteomes of marine microorganisms. ISME J. 6, 71–80 (2012).Margalef, R. Life-forms of phytoplankton as survival alternatives in an unstable environment. Ocean. Acta 1, (1978).Falkowski, P. G. & Oliver, M. J. Mix and match: How climate selects phytoplankton. Nat. Rev. Microbiol. 5, 813–819 (2007).Mende, D. R. et al. Environmental drivers of a microbial genomic transition zone in the ocean’s interior. Nat. Microbiol. 2, 1367–1373 (2017).Phoma, B. S. & Makhalanyane, T. P. Depth-dependent variables shape community structure and functionality in the Prince Edward Islands. Microb. Ecol. 81, 396–409 (2021).Kahru, M. & Mitchell, B. G. Seasonal and nonseasonal variability of satellite-derived chlorophyll and colored dissolved organic matter concentration in the California Current. J. Geophys. Res. Ocean. 106, 2517–2529 (2001).ADS 
CAS 
Article 

Google Scholar 
Barth, A., Walter, R. K., Robbins, I. & Pasulka, A. Seasonal and interannual variability of phytoplankton abundance and community composition on the Central Coast of California. Mar. Ecol. Prog. Ser. 637, (2020).Powell, J. R. & Ohman, M. D. Changes in zooplankton habitat, behavior, and acoustic scattering characteristics across glider-resolved fronts in the Southern California Current System. Prog. Oceanogr. 134, 77–92 (2015).ADS 
Article 

Google Scholar 
Taylor, A. G. & Landry, M. R. Phytoplankton biomass and size structure across trophic gradients in the southern California Current and adjacent ocean ecosystems. Mar. Ecol. Prog. Ser. 592, 1–17 (2018).ADS 
CAS 
Article 

Google Scholar 
Dutkiewicz, S., Follows, M. J. & Bragg, J. G. Modeling the coupling of ocean ecology and biogeochemistry. Glob. Biogeochem. Cycles 23, 1–15 (2009).Article 
CAS 

Google Scholar 
D’Ovidio, F., De Monte, S., Alvain, S., Dandonneau, Y. & Lévy, M. Fluid dynamical niches of phytoplankton types. Proc. Natl Acad. Sci. U. S. A. 107, 18366–18370 (2010).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Clayton, S., Dutkiewicz, S., Jahn, O. & Follows, M. J. Dispersal, eddies, and the diversity of marine phytoplankton. Limnol. Oceanogr. Fluids Environ. 3, 182–197 (2013).Article 

Google Scholar 
Moisan, T. A., Rufty, K. M., Moisan, J. R. & Linkswiler, M. A. Satellite observations of phytoplankton functional type spatial distributions, phenology, diversity, and ecotones. Front. Mar. Sci. 4, 1–24 (2017).Article 

Google Scholar 
Combes, V. et al. Cross-shore transport variability in the California Current: Ekman upwelling vs. eddy dynamics. Prog. Oceanogr. 109, 78–89 (2013).ADS 
Article 

Google Scholar 
Chenillat, F., Rivière, P., Capet, X., Franks, P. J. S. & Blanke, B. California coastal upwelling onset variability: cross-shore and bottom-up propagation in the planktonic ecosystem. PLoS ONE 8, (2013).Chenillat, F., Franks, P. J. S. & Combes, V. Biogeochemical properties of eddies in the California Current System. Geophys. Res. Lett. 43, 5812–5820 (2016).ADS 
CAS 
Article 

Google Scholar 
Edwards, K. F., Thomas, M. K., Klausmeier, C. A. & Litchman, E. Allometric scaling and taxonomic variation in nutrient utilization traits and maximum growth rate of phytoplankton. Limnol. Oceanogr. 57, 554–566 (2012).ADS 
Article 

Google Scholar 
Wells, B. K. et al. State of the California Current 2016–17: Still anything but ‘normal’ in the north. Calif. Cooperative Ocean. Fish. Investig. Rep. 58 (2017).Thompson, A. R. et al. State of the California Current 2017–18: Still not quite normal in the north and getting interesting in the south. Calif. Cooperative Ocean. Fish. Investig. Rep. 59 (2018).Ward, C. S. et al. Annual community patterns are driven by seasonal switching between closely related marine bacteria. ISME J. 11, 1412–1422 (2017).Bograd, S. J., Schroeder, I. D. & Jacox, M. G. A water mass history of the Southern California current system. Geophys. Res. Lett. 46, 6690–6698 (2019).ADS 
Article 

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 (2016).Amaral-Zettler, L. A., McCliment, E. A., Ducklow, H. W. & Huse, S. M. A method for studying protistan diversity using massively parallel sequencing of V9 hypervariable regions of small-subunit ribosomal RNA Genes. PLoS ONE 4, (2009).Bolyen, E. et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.J. 17, (2011).Callahan, B. J., Mcmurdie, P. J., Rosen, M. J., Han, A. W. & A, A. J. DADA2: High resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bokulich, N. A. et al. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome 6 (2018).Pedregosa, F. et al. Scikit-learn: Machine learning in Python. J. Mach. Learn. Res. 12 (2011).Pruesse, E. et al. SILVA: A comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids Res. 35 (2007).Guillou, L. et al. The Protist Ribosomal Reference database (PR2): A catalog of unicellular eukaryote Small Sub-Unit rRNA sequences with curated taxonomy. Nucleic Acids Res. 41 (2013).McMurdie, P. J. & Holmes, S. Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible. PLoS Comput. Biol. 10 (2014).Gloor, G. B., Wu, J. R., Pawlowsky-Glahn, V. & Egozcue, J. J. It’s all relative: analyzing microbiome data as compositions. Ann. Epidemiol. 26 (2016).Cameron, E. S., Schmidt, P. J., Tremblay, B. J. M., Emelko, M. B. & Müller, K. M. To rarefy or not to rarefy: Enhancing microbial community analysis through next-generation sequencing. bioRxiv. https://doi.org/10.1101/2020.09.09.290049 (2020).Oksanen, J. et al. vegan: Community Ecology Package. R package version 2.5-7. (2020).Bowman, J. S., Amaral-zettler, L. A., Rich, J. J., Luria, C. M. & Ducklow, H. W. Bacterial community segmentation facilitates the prediction of ecosystem function along the coast of the western Antarctic Peninsula. Nat. Publ. Gr. 11, 1460–1471 (2017).
Google Scholar 
Boelaert, J., Bendhaiba, L., Olteanu, M. & Villa-Vialaneix, N. SOMbrero: An R package for numeric and non-numeric self-organizing maps. Adv. Intell. Syst. Comput 295, 219–228 (2014).
Google Scholar 
Johnson, J. B. & Omland, K. S. Model selection in ecology and evolution. Trends Ecol. Evol. 19, 101–108 (2004).PubMed 
Article 

Google Scholar 
James, C. C. et al. Influence of nutrient supply on plankton microbiome biodiversity and distribution in a coastal upwelling region. https://doi.org/10.5281/zenodo.6359865 (2022).Legendre, P. & Legendre, L. Numerical ecology (Elsevier, 2012). More

Miyata, T., Miyazawa, S. & Yasunaga, T. Two types of amino acid substitutions in protein evolution. J. Mol. Evol. 12, 219–236 (1979).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Li, W.-H., Wu, C.-I. & Luo, C.-C. A new method for estimating synonymous and nonsynonymous rates of nucleotide substitution considering the relative likelihood of nucleotide and codon changes. Mol. Biol. Evol. 2, 150–174 (1985).PubMed 

Google Scholar 
Bielawski, J. P. & Yang, Z. Positive and negative selection in the DAZ gene family. Mol. Biol. Evol. 18, 523–529 (2001).CAS 
PubMed 
Article 

Google Scholar 
Ohta, T. Slightly deleterious mutant substitutions in evolution. Nature 246, 96–98 (1973).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Ohta, T. The nearly neutral theory of molecular evolution. Annu. Rev. Ecol. Evol. Syst. 23, 263–286 (1992).Article 

Google Scholar 
Johnson, K. P. & Seger, J. Elevated rates of nonsynonymous substitution in island birds. Mol. Biol. Evol. 18, 874–881 (2001).CAS 
PubMed 
Article 

Google Scholar 
Woolfit, M. & Bromham, L. Population size and molecular evolution on islands. Proc. Biol. Sci. 272, 2277–2282 (2005).CAS 
PubMed 
PubMed Central 

Google Scholar 
Ross, L., Hardy, N. B., Okusu, A. & Normark, B. B. Large population size predicts the distribution of asexuality in scale insects. Evolution 67, 196–206 (2013).PubMed 
Article 

Google Scholar 
Weber, C. C., Nabholz, B., Romiguier, J. & Ellegren, H. Kr/Kc but not dN/dS correlates positively with body mass in birds, raising implications for inferring lineage-specific selection. Genome Biol. 15, 542 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Brandt, A. et al. Effective purifying selection in ancient asexual oribatid mites. Nat. Commun. 8, 873 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Figuet, E. et al. Life history traits, protein evolution, and the nearly neutral theory in amniotes. Mol. Biol. Evol. 33(6), 1517–1527 (2016).CAS 
PubMed 
Article 

Google Scholar 
Saclier, N. et al. Life history traits impact the nuclear rate of substitution but not the mitochondrial rate in isopods. Mol. Biol. Evol. 35, 2900–2912 (2018).CAS 
PubMed 
Article 

Google Scholar 
Hebert, P. D. The Daphnia of North America: An Illustrated Fauna (on CD-ROM) (CyberNatural Software, Guelph, 1995).
Google Scholar 
Colbourne, J. K. et al. Phylogenetics and evolution of a circumarctic species complex (Cladocera: Daphnia pulex). Biol. J. Linn. Soc. 65, 347–365 (1998).
Google Scholar 
Crease, T. J., Omilian, A. R., Costanzo, K. S. & Taylor, D. J. Transcontinental phylogeography of the Daphnia pulex species complex. PLoS ONE 7, e46620 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Mergeay, J., Verschuren, D. & De Meester, L. Cryptic invasion and dispersal of an American Daphnia in East Africa. Limnol. Oceanogr. 50, 1278–1283 (2005).ADS 
CAS 
Article 

Google Scholar 
Ma, X. et al. Lineage diversity and reproductive modes of the Daphnia pulex group in Chinese lakes and reservoirs. Mol. Phylogenet. Evol. 130, 424–433 (2019).PubMed 
Article 

Google Scholar 
So, M. et al. Invasion and molecular evolution of Daphnia pulex in Japan. Limnol. Oceanogr. 60, 1129–1138 (2015).ADS 
Article 

Google Scholar 
Duggan, I. C. et al. Identifying invertebrate invasions using morphological and molecular analyses: North American Daphnia ‘pulex’ in New Zealand fresh waters. Aquat. Invasions 7, 585–590 (2012).Article 

Google Scholar 
Ye, Z. et al. The rapid, mass invasion of New Zealand by North American Daphnia “pulex”. Limnol. Oceanogr. 66, 2673–2683 (2021).ADS 
Article 

Google Scholar 
Paland, S., Colbourne, J. K. & Lynch, M. Evolutionary history of contagious asexuality in Daphnia pulex. Evolution 59, 800–813 (2005).CAS 
PubMed 
Article 

Google Scholar 
Muller, H. J. The relation of recombination to mutational advance. Mutat. Res. 106, 2–9 (1964).CAS 
PubMed 
Article 

Google Scholar 
Felsenstein, J. The evolutionary advantage of recombination. Genetics 78, 737–756 (1974).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Paland, S. & Lynch, M. Transitions to asexuality result in excess amino acid substitutions. Science 311, 990–992 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Johnson, S. G. & Howard, R. S. Contrasting patterns of synonymous and nonsynonymous sequence evolution in asexual and sexual freshwater snail lineages. Evolution 61, 2728–2735 (2007).PubMed 
Article 

Google Scholar 
Neiman, M. et al. Accelerated mutation accumulation in asexual lineages of a freshwater snail. Mol. Biol. Evol. 27, 954–963 (2010).CAS 
PubMed 
Article 

Google Scholar 
Henry, L., Schwander, T. & Crespi, B. J. Deleterious mutation accumulation in asexual Timema stick insects. Mol. Biol. Evol. 29, 401–408 (2012).CAS 
PubMed 
Article 

Google Scholar 
Tucker, A. E. et al. Population-genomic insights into the evolutionary origin and fate of obligately asexual Daphnia pulex. Proc. Natl. Acad. Sci. 110, 15740–15745 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Colbourne, J. K. et al. The ecoresponsive genome of Daphnia pulex. Science 331, 555–561 (2011).ADS 
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 
Ye, Z. et al. A new reference genome assembly for the microcrustacean Daphnia pulex. G3 (Bethesda) 7, 1405–1416 (2017).CAS 
Article 

Google Scholar 
Keith, N. et al. High mutational rates of large-scale duplication and deletion in Daphnia pulex. Genome Res. 26, 60–69 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hall, D. J. An experimental approach to the dynamics of a natural population of Daphnia galeata mendotae. Ecology 45, 94–112 (1964).Article 

Google Scholar 
McCauley, E., Murdoch, W. W. & Nisbet, R. M. Growth, reproduction, and mortality of Daphnia pulex Leydig: Life at low food. Ecology 4, 505–514 (1990).
Google Scholar 
Xu, S. et al. High mutation rates in the mitochondrial genomes of Daphnia pulex. Mol. Biol. Evol. 29, 763–769 (2012).CAS 
PubMed 
Article 

Google Scholar 
Zheng, Y., Peng, R., Kuro-o, M. & Zeng, X. Exploring patterns and extent of bias in estimating divergence time from mitochondrial DNA sequence data in a particular lineage: A case study of salamanders (Order Caudata). Mol. Biol. Evol. 28, 2521–2535 (2011).CAS 
PubMed 
Article 

Google Scholar 
Zaret, T. M. Predation and Freshwater Communities (Yale University Press, New Haven, 1980).
Google Scholar 
Lynch, M. Predation, competition, and zooplankton community structure: An experimental study. Limnol. Oceanogr. 24, 253–272 (1979).ADS 
Article 

Google Scholar 
Mills, E. L. & Forney, J. L. Impact on Daphnia pulex of predation by young yellow perch in Oneida Lake, New York. Trans. Am. Fish. Soc. 112(2A), 154–161 (1983).Article 

Google Scholar 
Craddock, D. R. Effects of increased water temperature on Daphnia pulex. Fish. Bull. 74, 403–408 (1976).
Google Scholar 
Maruoka, N. & Urabe, J. Inter and intraspecific competitive abilities and the distribution ranges of two Daphnia species in Eurasian continental islands. Popul. Ecol. 62, 353–363 (2020).Article 

Google Scholar 
Dodson, S. I. & Hanazato, T. Commentary on effects of anthropogenic and natural organic chemicals on development, swimming behavior, and reproduction of Daphnia, a key member of aquatic ecosystems. Environ. Health Perspect. 103(Suppl 4), 7–11 (1995).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Claska, M. E. & Gilbert, J. J. The effect of temperature on the response of Daphnia to toxic cyanobacteria. Freshw. Biol. 39, 221–232 (1998).Article 

Google Scholar 
Bast, J. et al. Consequences of asexuality in natural populations: Insights from stick insects. Mol. Biol. Evol. 35, 1668–1677 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hartfield, M. Evolutionary genetic consequences of facultative sex and outcrossing. J Evol Biol 29, 5–22 (2016).CAS 
PubMed 
Article 

Google Scholar 
Hörandl, E. et al. Genome evolution of asexual organisms and the paradox of sex in eukaryotes. In Evolutionary Biology—A Transdisciplinary Approach (ed. Pontarotti, P.) (Springer, Cham, 2020). https://doi.org/10.1007/978-3-030-57246-4_7.Chapter 

Google Scholar 
Lynch, M., Bürger, R., Butcher, D. & Gabriel, W. The mutational meltdown in asexual populations. J. Hered. 84, 339–344 (1993).CAS 
PubMed 
Article 

Google Scholar 
Gordo, I. & Charlesworth, B. The degeneration of asexual haploid populations and the speed of Muller’s ratchet. Genetics 154, 1379–1387 (2000).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Downing, J. A. et al. The global abundance and size distribution of lakes, ponds, and impoundments. Limnol. Oceanogr. 51, 2388–2397 (2006).ADS 
Article 

Google Scholar 
McDonald, C. P., Rover, J. A., Stets, E. G. & Striegl, R. G. The regional abundance and size distribution of lakes and reservoirs in the United States and implications for estimates of global lake extent. Limnol. Oceanogr. 57, 597–606 (2012).ADS 
Article 

Google Scholar 
De Meester, L., Góme, A., Okamura, B. & Schwenk, K. The monopolization hypothesis and the dispersal-gene flow paradox in aquatic organisms. Acta Oecol. 23, 121–135 (2002).ADS 
Article 

Google Scholar 
Fukami, T., Bezemer, T. M., Mortimer, S. R. & Van Der Putten, W. H. Species divergence and trait convergence in experimental plant community assembly. Ecol. Lett. 8, 1283–1290 (2005).Article 

Google Scholar 
Makino, T. & Kawata, M. Invasive invertebrates associated with highly duplicated gene content. Mol. Ecol. 28, 1652–1663 (2019).CAS 
PubMed 
Article 

Google Scholar 
Kondrashov, F. A. Gene duplication as a mechanism of genomic adaptation to a changing environment. Proc. R. Soc. Lond. B Biol. Sci. 279, 5048–5057 (2012).
Google Scholar 
Barrick, J. E. & Lenski, R. E. Genome dynamics during experimental evolution. Nat. Rev. Genet. 14, 827–839 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Rocha, E. P. C. Neutral theory, microbial practice: Challenges in bacterial population genetics. Mol. Biol. Evol. 35, 1338–1347 (2018).CAS 
PubMed 
Article 

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 
Article 

Google Scholar 
Tanabe, A. S. Kakusan4 and Aminosan: Two programs for comparing nonpartitioned, proportional and separate models for combined molecular phylogenetic analyses of multilocus sequence data. Mol. Ecol. Resour. 11, 914–921 (2011).PubMed 
Article 

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

Google Scholar 
Tian, X., Ohtsuki, H. & Urabe, J. Evolution of asexual Daphnia pulex in Japan: Variations and covariations of the digestive, morphological and life history traits. BMC Evol. Biol. 19, 122 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Chen, Y. et al. SOAPnuke: A MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience 7, 1–6 (2018).ADS 
PubMed 
PubMed Central 

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

Google Scholar 
Lee, T. H. et al. SNPhylo: A pipeline to construct a phylogenetic tree from huge SNP data. BMC Genomics 15, 162 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
R Core Team, R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, Austria, 2019). https://www.R-project.org/ More

Meropol, S. B., Haupt, A. A. & Debanne, S. M. Incidence and outcomes of infections caused by multidrug-resistant Enterobacteriaceae in Children, 2007–2015. J. Pediatr. Infect. Dis. Soc. 7, 36–45 (2018).Article 

Google Scholar 
Moxon, C. A. & Paulus, S. Beta-lactamases in Enterobacteriaceae infections in children. J. Infect. 72, S41–S49 (2016).PubMed 
Article 

Google Scholar 
Morrissey, I. et al. A review of ten years of the study for monitoring antimicrobial resistance trends (SMART) from 2002 to 2011. Pharmaceuticals 6, 1335–1346 (2013).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Junnila, J. et al. Changing epidemiology of methicillin-resistant Staphylococcus aureus in a low endemicity area—new challenges for MRSA control. Eur. J. Clin. Microbiol. Infect. Dis. 39, 2299–2307 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Milstone, A. M. et al. Methicillin-resistant Staphylococcus aureus colonization and risk of subsequent infection in critically ill children: Importance of preventing nosocomial methicillin-resistant Staphylococcus aureus transmission. Clin. Infect. Dis. 53, 853–859 (2011).PubMed 
PubMed Central 
Article 

Google Scholar 
Lakhundi, S. & Zhang, K. Methicillin-resistant Staphylococcus aureus: Molecular characterization, evolution, and epidemiology. Clin. Microbiol. Rev. 31, e00020-18 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Schlesinger, Y. et al. Methicillin-resistant nasal colonization in children in Jerusalem: Community vs. chronic care institutions. Isr. Med. Assoc. J. 5, 847–851 (2003).PubMed 

Google Scholar 
Liang, B. et al. Active surveillance, drug resistance, and genotypic profiling of Staphylococcus aureus among school-age children in China. Front. Med. 8, 701494 (2021).Article 

Google Scholar 
Del Rosal, T. et al. Staphylococcus aureus nasal colonization in Spanish children. The COSACO Nationwide Surveillance Study. Infect. Drug Resist. 13, 4643–4651 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Feudtner, C., Feinstein, J. A., Zhong, W., Hall, M. & Dai, D. Pediatric complex chronic conditions classification system version 2: Updated for ICD-10 and complex medical technology dependence and transplantation. BMC Pediatr. 14, 199 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Climent Alcalá, F. J., García Fernández de Villalta, M., Escosa García, L., Rodríguez Alonso, A. & Albajara Velasco, L. A. Unidad de niños con patología crónica compleja. Un modelo necesario en nuestros hospitales. Anales de Pediatría 88, 12–18 (2018).PubMed 
Article 

Google Scholar 
Gesualdo, F. et al. Methicillin-resistant Staphylococcus aureus nasal colonization in a department of pediatrics: A cross-sectional study. Ital. J. Pediatr. 40, 3 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Yamamoto, M. et al. Effective surveillance to identify the surgical patients carrying methicillin-resistant Staphylococcus aureus on admission in a pediatric ward. Osaka City Med. J. 62, 1–9 (2016).PubMed 

Google Scholar 
Lukac, P. J., Bonomo, R. A. & Logan, L. K. Extended-spectrum-lactamase-producing Enterobacteriaceae in children: Old foe, emerging threat. Clin. Infect. Dis. https://doi.org/10.1093/cid/civ020 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Fedler, K. A., Biedenbach, D. J. & Jones, R. N. Assessment of pathogen frequency and resistance patterns among pediatric patient isolates: Report from the 2004 SENTRY Antimicrobial Surveillance Program on 3 continents. Diagn. Microbiol. Infect. Dis. 56, 427–436 (2006).CAS 
PubMed 
Article 

Google Scholar 
Caselli, D. et al. Incidence of colonization and bloodstream infection with carbapenem-resistant Enterobacteriaceae in children receiving antineoplastic chemotherapy in Italy. Infect. Dis. 48, 152–155 (2016).Article 

Google Scholar 
Logan, L. K. et al. Multidrug- and Carbapenem-Resistant Pseudomonas aeruginosa in Children, United States, 1999–2012. JPIDSJ piw064 (2016) https://doi.org/10.1093/jpids/piw064.Flokas, M. E., Alevizakos, M., Shehadeh, F., Andreatos, N. & Mylonakis, E. Extended-spectrum β-lactamase-producing Enterobacteriaceae colonisation in long-term care facilities: A systematic review and meta-analysis. Int. J. Antimicrob. Agents 50, 649–656 (2017).CAS 
PubMed 
Article 

Google Scholar 
Bharadwaj, R. et al. Drug-resistant Enterobacteriaceae colonization is associated with healthcare utilization and antimicrobial use among inpatients in Pune, India. BMC Infect. Dis. 18, 504 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Logan, L. K. Carbapenem-resistant Enterobacteriaceae: An emerging problem in children. Clin. Infect. Dis. 55, 852–859 (2012).CAS 
PubMed 
Article 

Google Scholar 
Badal, R. E. et al. Etiology, extended-spectrum β-lactamase rates and antimicrobial susceptibility of gram-negative bacilli causing intra-abdominal infections in patients in general pediatric and pediatric intensive care units—global data from the Study for Monitoring Antimicrobial Resistance Trends 2008 to 2010. Pediatr. Infect. Dis. J. 32, 636–640 (2013).PubMed 
Article 

Google Scholar 
Wang, Q. et al. Risk factors and clinical outcomes for carbapenem-resistant Enterobacteriaceae nosocomial infections. Eur. J. Clin. Microbiol. Infect. Dis. 35, 1679–1689 (2016).CAS 
PubMed 
Article 

Google Scholar 
Sahbudak Bal, Z. et al. The prospective evaluation of risk factors and clinical influence of carbapenem resistance in children with gram-negative bacteria infection. Am. J. Infect. Control 46, 147–153 (2018).CAS 
PubMed 
Article 

Google Scholar 
Simon, T. D. et al. Pediatric medical complexity algorithm: A new method to stratify children by medical complexity. Pediatrics 133, e1647–e1654 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Román, F. et al. Characterization of methicillin-resistant Staphylococcus aureus strains colonizing the nostrils of Spanish children. MicrobiologyOpen 10, e1235 (2021).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
EUCAST. European committee on antimicrobial susceptibility testing breakpoint tables for interpretation of MICs and zone diameters. The European Committee on Antimicrobial Susceptibility Testing. (2018).Oteo, J. et al. Prospective multicenter study of carbapenemase-producing Enterobacteriaceae from 83 hospitals in Spain reveals high in vitro susceptibility to colistin and meropenem. Antimicrob. Agents Chemother. 59, 3406–3412 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Maseda, E. et al. Risk factors for colonization by carbapenemase-producing enterobacteria at admission to a Surgical ICU: A retrospective study. Enferm. Infecc. Microbiol. Clin. 35, 333–337 (2017).PubMed 
Article 

Google Scholar 
Bassetti, M., Nicco, E. & Mikulska, M. Why is community-associated MRSA spreading across the world and how will it change clinical practice?. Int. J. Antimicrob. Agents 34, S15–S19 (2009).CAS 
PubMed 
Article 

Google Scholar 
El Cheikh, M. R., Barbosa, J. M., Caixêta, J. A. S. & Avelino, M. A. G. Microbiology of tracheal secretions: What to expect with children and adolescents with tracheostomies. Int. Arch. Otorhinolaryngol. 22, 50–54 (2018).PubMed 
Article 

Google Scholar 
González-Del Castillo, J. et al. BAHNG score: Predictive model for detection of subjects with the oropharynx colonized by uncommon microorganisms. Rev. Esp Quimioter. 30, 422–428 (2017).PubMed 

Google Scholar 
Hu, X. et al. Risk factors for methicillin-resistant Staphylococcus aureus colonization and infection in patients with human immunodeficiency virus infection: A systematic review and meta-analysis. J. Int. Med. Res. 50, 3000605211063019 (2022).CAS 
PubMed 

Google Scholar 
Gleeson, A., Larkin, P., Walsh, C. & O’Sullivan, N. Methicillin-resistant Staphylococcus aureus: Prevalence, incidence, risk factors, and effects on survival of patients in a specialist palliative care unit: A prospective observational study. Palliat. Med. 30, 374–381 (2016).PubMed 
Article 

Google Scholar 
Hogardt, M. et al. Current prevalence of multidrug-resistant organisms in long-term care facilities in the Rhine-Main district, Germany, 2013. Euro Surveill. 20, 21171 (2015).PubMed 
Article 

Google Scholar 
Warren, D. K. et al. Epidemiology of methicillin-resistant Staphylococcus aureus colonization in a surgical intensive care unit. Infect. Control Hosp. Epidemiol. 27, 1032–1040 (2006).PubMed 
Article 

Google Scholar 
Folgori, L. et al. Healthcare-associated infections in pediatric and neonatal intensive care units: Impact of underlying risk factors and antimicrobial resistance on 30-day case-fatality in Italy and Brazil. Infect. Control Hosp. Epidemiol. 37, 1302–1309 (2016).PubMed 
Article 

Google Scholar 
Béranger, A. et al. Early bacterial infections after pediatric liver transplantation in the era of multidrug-resistant bacteria: Nine-year single-center retrospective experience. Pediatr. Infect. Dis. J. 39, e169–e175 (2020).PubMed 
Article 

Google Scholar 
Bouras, D. et al. Staphylococcus aureus osteoarticular infections in children: An 8-year review of molecular microbiology, antibiotic resistance and clinical characteristics. J. Med. Microbiol. 67, 1753–1760 (2018).MathSciNet 
CAS 
PubMed 
Article 

Google Scholar 
Rodriguez, M., Hogan, P. G., Krauss, M., Warren, D. K. & Fritz, S. A. Measurement and impact of Staphylococcus aureus colonization pressure in households. J. Pediatr. Infect. Dis. Soc. 2, 147–154 (2013).Article 

Google Scholar 
Messina, N. L., Williamson, D. A., Robins-Browne, R., Bryant, P. A. & Curtis, N. Risk factors for carriage of antibiotic-resistant bacteria in healthy children in the community: A systematic review. Pediatr. Infect. Dis. J. 39, 397–405 (2020).PubMed 
Article 

Google Scholar 
Dualleh, N. et al. Colonization with multiresistant bacteria in acute hospital care: The association of prior antibiotic consumption as a risk factor. J. Antimicrob. Chemother. 75, 3675–3681 (2020).CAS 
PubMed 
Article 

Google Scholar 
Daskalaki, M. et al. Panton-Valentine leukocidin-positive Staphylococcus aureus skin and soft tissue infections among children in an emergency department in Madrid, Spain. Clin. Microbiol. Infect. 16, 74–77 (2010).CAS 
PubMed 
Article 

Google Scholar 
Aguilera-Alonso, D., Escosa-García, L., Saavedra-Lozano, J., Cercenado, E. & Baquero-Artigao, F. Carbapenem-resistant gram-negative bacterial infections in children. Antimicrob. Agents Chemother. 64, e02183-e2219 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Phichaphop, C. et al. High prevalence of multidrug-resistant gram-negative bacterial infection following pediatric liver transplantation. Medicine 99, e23169 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Tacconelli, E. et al. ESCMID guidelines for the management of the infection control measures to reduce transmission of multidrug-resistant Gram-negative bacteria in hospitalized patients. Clin. Microbiol. Infect. 20, 1–55 (2014).PubMed 
Article 

Google Scholar 
McConville, T. H., Sullivan, S. B., Gomez-Simmonds, A., Whittier, S. & Uhlemann, A.-C. Carbapenem-resistant Enterobacteriaceae colonization (CRE) and subsequent risk of infection and 90-day mortality in critically ill patients, an observational study. PLoS ONE 12, e0186195 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Tamma, P. D. et al. The likelihood of developing a carbapenem-resistant Enterobacteriaceae Infection during a hospital stay. Antimicrob. Agents Chemother. 63, e00757-e819 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Detsis, M., Karanika, S. & Mylonakis, E. ICU acquisition rate, risk factors, and clinical significance of digestive tract colonization with extended-spectrum beta-lactamase-producing Enterobacteriaceae: A systematic review and meta-analysis. Crit. Care Med. 45, 705–714 (2017).PubMed 
Article 

Google Scholar  More

The haplotype networks, PCoA analysis and STRUCTURE analysis based on the six nuclear loci confirm that S. longifolium is a hybrid between S. emersum and S. gramineum, providing molecular support for previous morphological analyses5. Furthermore, all individuals with intermediate admixture coefficient (Fig. 2b) and private haplotypes only present in one out of six nuclear loci (Fig. 1) suggest that S. longifolium is most likely a F1 hybrid. We thus hypothesized that S. emersum and S. gramineum could likely maintain their species boundary through the post-zygote reproductive isolation mechanism of F1 generation sterility. This hypothesis is possible based on the observations from hybrids in European Russia. The pollen viability was checked in S. longifolium samples from Vysokovskoe Lake and Sabro Lake, and the vast majority of checked pollens were sterile5. In addition, flowering plants of S. longifolium often do not form seeds, or the seeds are puny and significantly inferior to normal seeds in size5. However, the hypothesis is only based on our limited sampling, which is contrary to the conclusion inferred from morphological characteristics that it is fertile and may backcross with parental species1. Further studies with extensive sampling are necessary to test our hypothesis.The chloroplast DNA fragment trnH-psbA was used to infer the direction of hybridization between S. emersum and S. gramineum because chloroplast DNA is maternal inheritance in Sparganium3,4. The hybrid S. longifolium shared haplotypes with S. emersum and S. gramineum simultaneously (Fig. 1). This finding clearly indicates that bidirectional hybridization exists between S. emersum and S. gramineum. At the same time, the different frequency of these two haplotypes in the hybrid (H1, 19.1% vs. H2, 80.9%) means that the direction of hybridization is asymmetric. A variety of factors can lead to asymmetry in natural hybridization, such as flowering time, preference of pollinators, quality and quantity of pollen, cross incompatibility and the abundance of parent species7,8. Rare species usually act as maternal species relative to abundant species9,10. S. gramineum is confined to oligotrophic lakes and its abundance is obviously lower than that of S. emersum1,11. The relatively scarcity combined with the ecology of S. gramineum make it more often act as maternal species when hybridizing with S. emersum.As described by5, the morphological diversification of S. longifolium was also observed in this study. For example, individuals of S. longifolium with emergent and floating-leaved life forms occur concurrently in Zaozer’ye Lake (Supplementary Fig. S2). However, all individuals had the same haplotype H2 as S. gramineum (Fig. 1), suggesting that the direction of hybridization do not determine life form of S. longifolium. In addition, all individuals of S. longifolium sampled here are likely F1 hybrid. Their variable phenotypes could not be associated with traits segregation due to F2 generation or backcross. Detailed ecological investigation combining with research at the genomic level are essential to find out the potential factors leading to morphological diversification of S. longifolium.Here, using sequences of six nuclear loci and one chloroplast DNA fragment, we confirmed that S. longifolium is the hybrid between S. emersum and S. gramineum. The natural hybridization between S. emersum and S. gramineum is bidirectional but the latter mainly acts as maternal species. We also found that all samples of S. longifolium were F1 generations, indicating that S. emersum and S. gramineum could maintain their species boundary through the post-zygote reproductive isolation mechanism of F1 generation sterility. More

Overall strategyAn admixture analysis aims to estimate the admixture proportions (or ancestries), Q, of each sampled individual in a given number of K source populations (Pritchard et al. 2000), and the characteristic allele frequencies, P, at each locus of each inferred source population. Even though Q is frequently of the primary interest, P must be estimated simultaneously because we have genotype data only and Q is highly dependent on P which actually defines the source populations. For N individuals from K source populations genotyped at L loci with a total number of A alleles, the numbers of independent variables in Q and P are VQ = (K − 1)N and VP = (A − L)K, respectively. The high dimensionality of an admixture analysis, with V = VQ + VP = (K − 1)N + (A − L)K variables, not only incurs a large computational burden, but also poses a high risk of non-convergence (to the global maximum) for any algorithm, especially when either Q or P is expected to be poorly estimated in difficult situations such as a small sample (say, a couple) of individuals from each source population or low differentiation.I propose a two-step procedure with corresponding algorithms to reduce the risk of non-convergence, to speed up the computation, and to make more accurate inferences of both Q and P. In the first step, I assume a mixture model (Pritchard et al. 2000; Falush et al. 2003) that individuals in a sample can come from different source populations, but each individual’s genome comes exclusively from a single population. Under this simplified probabilistic model, I conduct a clustering analysis to obtain estimates of both individual memberships and allele frequencies of each cluster by a global maximisation algorithm, simulated annealing, with extra care (details below) of convergence. In the absence of admixture and with sufficient information for complete recovery of population structure, the estimated individual memberships and allele frequencies of the clusters are expected to be equivalent to Q (with element qik = 1 and qil = 0 if individual i is inferred to be in cluster k where l ≠ k) and P, respectively. Otherwise, they are expected to be good approximations of Q and P, because an admixed individual i with the highest ancestral proportion from a population would be expected to be assigned (exclusively) to that population. In the second step, I assume an admixture model (Pritchard et al. 2000; Falush et al. 2003) to refine estimates of Q and P, using an EM algorithm and the start parameter (Q and P) values obtained from the clustering analysis. Because the starting values are already close to the truth, the algorithm is fast and has a much-reduced risk of converging to a local maximum than the original EM algorithms (Tang et al. 2005; Alexander et al. 2009).Clustering analysisI assume N diploid individuals are sampled from K source populations. The origin of a sampled individual from the K source populations is unknown, which is the primary interest of structure analysis. However, if it is (partially) known, this information can be used to supervise (help) the clustering analysis of other sampled individuals of unknown origins. Each individual’s genome comes exclusively from one of the K unknown source populations (i.e., mixture model, no admixture). I assume each individual is genotyped at L loci, with a diploid genotype {xil1, xil2} for individual i (=1, 2, …, N) at locus l (=1, 2, …, L). The task of the clustering analysis is to sort the N individuals with genotype data X = {xila:i = 1, 2, …, N; l = 1, 2, …, L; a = 1, 2} into K clusters, with each representing a source population. No assumption is made about the evolutionary relationships of the populations, which, when summarized by F statistics, are estimated from the same genotype data in both clustering and admixture analyses.Suppose, in a given clustering configuration Ω = {Ω1, Ω2, …, ΩK}, cluster k (=1, 2, …, K), Ωk, contains a set of Nk (with Nk  > 0 and (mathop {sum}nolimits_{k = 1}^K {N_k equiv N})) individuals, denoted by Ωk = {ωk1, ωk2, …, ωkNk} where ωkj is the index of the jth individual in cluster k. The genotype data of the Nk individuals in cluster k is Xk = {xila: i ∈ Ωk; l = 1, 2, …, L; a = 1, 2}. The log-likelihood of Ωk is then the log probability of observing Xk given Ωk$${{{mathcal{L}}}}_kleft( {{{{mathbf{Omega }}}}_k} right) = {{{mathrm{LogP}}}}left( {{{{mathbf{X}}}}_kleft| {{{{mathbf{Omega }}}}_k} right.} right) = mathop {sum}limits_{l = 1}^L {mathop {sum}limits_{j = 1}^{J_l} {c_{klj}{{{mathrm{Log}}}}left( {p_{klj}} right)} }$$
(1)
where cklj and pklj are the count of copies and the frequency, respectively, of allele j at locus l in cluster k, and Jl is the number of alleles at locus l. Given Ωk, cklj is counted from genotype data Xk, and allele frequency pklj is estimated by$$p_{klj} = left( {p_{lj} + c_{klj}} right)/mathop {sum}limits_{m = 1}^{J_l} {left( {p_{lm} + c_{klm}} right)}$$
(2)
where plj is the frequency of allele j at locus l in the entire population represented by the K clusters. plj is calculated by$$p_{lj} = mathop {sum}limits_{k = 1}^K {c_{klj}} /mathop {sum}limits_{m = 1}^{J_l} {mathop {sum}limits_{k = 1}^K {c_{klm}} } = c_{lj}/mathop {sum}limits_{m = 1}^{J_l} {c_{lm}}$$
(3)
where (c_{lm} = mathop {sum}nolimits_{k = 1}^K {c_{klm}}) is the count of allele m (=1, 2, …, Jl) at locus l in the entire sample of individuals.Under the mixture model above, clusters are only weakly dependent (with the extent of dependency decreasing with an increasing value of K) and the total log-likelihood of the clustering configuration, Ω = {Ω1, Ω2, …, ΩK}, is thus$${{{mathcal{L}}}}left( {{{mathbf{Omega }}}} right) = mathop {sum}limits_{k = 1}^K {{{{mathcal{L}}}}_kleft( {{{{mathbf{Omega }}}}_k} right)} ,$$
(4)
where ({{{mathcal{L}}}}_kleft( {{{{mathbf{Omega }}}}_k} right)) is calculated by (1).It is worth noting that allele frequencies, P, are modelled as hidden or nuisance variables and are estimated as a by-product of maximising (4) for estimates of Ω. Yet, careful modelling of P proves important for estimating Ω, as the two are highly dependent. Bayesian admixture methods assume allele frequencies pkl = {pkl1, pkl2, …, (p_{klj_l})} in a Dirichlet distribution (e.g., Foreman et al. 1997; Rannala and Mountain 1997; Pritchard et al. 2000), ({{{mathcal{D}}}}left( {lambda _1,lambda _2, ldots ,lambda _{J_l}} right)). For any population k, the uncorrelated (Pritchard et al. 2000) and correlated (Falush et al. 2003) allele frequency model assumes λj = 1 and (lambda_j=p_{ol_j}F_K/(1-F_k)), respectively, for j = 1, 2, …, Jl. In the latter model, p0lj is the frequency of allele j at locus l in the ancestral population (common to the K derived populations), and Fk is the differentiation of population k from the ancestral population. In contrast, likelihood admixture methods (e.g., Tang et al. 2005; Alexander et al. 2009; Frichot et al. 2014) and non-model based clustering methods (e.g., K-means method, Jombart et al. 2010) do not use any prior, which is equivalent to assuming plj ≡ 0 for j = 1, 2, …, Jl in Eq. (2). However, properly modelling prior allele frequencies, as carefully considered in Bayesian methods (Pritchard et al. 2000; Falush et al. 2003), becomes important in situations where allele frequencies are not well defined or tricky to estimate, such as when few individuals are sampled from a source population or when rare alleles are present. The frequentist estimator (2) is in spirit similar to the Bayesian correlated allele frequency model (Falush et al. 2003), and leads to accurate results in various situations to be shown in this study. I have also tried alternatives such as plj ≡ 1/Jl (which is similar to the uncorrelated allele frequency model of Pritchard et al. 2000) or plj ≡ 0 (which is equivalent to the treatment in previous likelihood admixture analysis or non-model based clustering analysis) in replacement of (2), but none works as well as (2) and could yield much less accurate results in difficult situations (below).Scaling for unbalanced samplingBayesian methods of STRUCTURE’s admixture model assume an individual i’s ancestry, qi = {qi1, qi2, …, qiK}, follows a prior Dirichlet probability distribution ({{{mathbf{q}}}}_isim {{{mathbf{{{{mathcal{D}}}}}}}}left( {alpha _1,alpha _2, ldots ,alpha _K} right)) (Pritchard et al. 2000; Falush et al. 2003). By default, α1 = α2 = ··· = αK = α, which essentially assumes that an individual has its ancestry originating from each of the assumed K populations at an equal prior probability of 1/K. To model unequal sample sizes such that an individual comes from a more intensively sampled population at a higher prior probability, STRUCTURE also has applied an alternative prior, α1 ≠ α2 ≠ ··· ≠ αK. It is shown that, when sampling intensity is heavily unbalanced among populations, the default prior could lead to the split of a large cluster and the merge of small clusters, while the alternative prior yields much more accurate results (Wang 2017). These priors have a large impact on admixture analysis; applying the default prior to data of highly unbalanced samples leads to inaccurate Q estimates even when many informative markers are used (Wang 2017).Unfortunately, current non-model based or likelihood-based admixture analysis methods do not utilise these or other priors for handling unbalanced sampling. As a result, they can give inaccurate admixture estimates, just like STRUCTURE under the default ancestry prior model, for data from highly unbalanced sampling. To reduce the cluster split and merge problems, herein I propose the following method to scale the likelihood of a cluster by the size, the number of individual members, of the cluster.The original log-likelihood of cluster k, ({{{mathcal{L}}}}_kleft( {{{{mathbf{Omega }}}}_k} right)), is calculated by (1). It is then scaled by the cluster size, Nk, as$${{{mathcal{L}}}}_{Sk}left( {{{{mathbf{Omega }}}}_k} right) = {{{mathcal{L}}}}_kleft( {{{{mathbf{Omega }}}}_k} right)/left( {1 + e^{sN_k/left( {8N} right)}} right),$$
(5)
where s is the scaling factor taking values 1, 2, 3 for weak, medium and strong scaling, respectively. This scaling scheme encourages large clusters and discourages small clusters. Although (5) is not an analytically derived but an empirical equation and is thus not guaranteed to be optimal, extensive simulations (some shown below) verify that the scaling scheme works very well for data from highly unbalanced sampling, yielding accurate clustering analysis results and thus similarly or more accurate admixture estimates than STRUCTURE under its alternative ancestry model. The most appropriate scaling level (1, 2 or 3) for a particular dataset depends on how unbalanced the sampling is, how much differentiated the populations are, and how much informative the markers are. For example, a low scaling level, s = 1, is appropriate when many markers are genotyped for a set of lowly differentiated (low FST) populations. Usually, we do not know these factors in analysing the data. Therefore, when the data are suspected to be unbalanced in sampling among populations, they are better analysed comparatively with different levels of scaling (0, 1, 2, and 3). When the applied level of scaling is too low, large populations tend to be split and small populations tend to be merged. When the applied level of scaling is too high, small populations tend to be merged among themselves or with a large population. With the help of some internal information such as consistency of replicate runs at the same scaling level and the same K value and some external information such as sampling locations in examining the admixture estimates, the appropriate scaling level can be determined.Simulated annealing algorithmA likelihood function with many variables, such as (4), is difficult to maximise for estimates of the variables. Traditional methods, such as derivative based Newton-Raphson algorithm (e.g., Tang et al. 2005) and non-derivative based EM algorithm (Dempster et al. 1977; Tang et al. 2005; Alexander et al. 2009), may converge to a local rather than the global maximum for a large scale problem with ridges and plateaus (Gaffe et al. 1994). Although trying multiple replicate runs with different starting values and choosing the run with the highest likelihood could reduce the risk of landing on a local maximum, a global maximum cannot be guaranteed regardless of the number of runs. The Bayesian approach as implemented in STRUCTURE (Pritchard et al. 2000) has a similar problem, as different replicate runs of the same data with the same parameter and model choices but different random number seeds may yield different admixture estimates and likelihood values (Tang et al. 2005; below).Simulated annealing (SA) was developed to optimise very large and complex systems (Kirkpatrick et al. 1983). Using the Metropolis algorithm (Metropolis et al. 1953) from statistical mechanics, SA can find the global maximum by searching both downhill and uphill and by traversing deep valleys on the likelihood surface to avoid getting stuck on a local maximum (Kirkpatrick et al. 1983; Goffe et al. 1994). It has been proved to be highly powerful in pedigree reconstruction (Wang 2004; Wang and Santure 2009) from genotype data, which is probably more difficult than population structure reconstruction (i.e., clustering analysis) because the genetic structure (i.e., sibship) of the former is, in general, more numerous, more complicated with hierarchy, and smaller (thus more elusive and more difficult to define) than that in the latter. Herein I propose a SA algorithm for a population clustering analysis, as detailed in Supplementary Appendix 1.Admixture analysisUnder the mixture model, the above clustering analysis partitions the N sampled individuals into a predefined K clusters, each representing a source population. The properties (e.g., genetic diversity) of and the relationships (e.g., FST) among these populations can be learnt from the inferred clusters. However, the clustering results are accurate only when the mixture model is valid. For a sample containing a substantial proportion of highly admixed individuals (i.e., who have recent ancestors from multiple source populations), the clustering results are just approximations. In such a case, the admixture model is more appropriate and can be used to refine the mixture analysis results by inferring the admixture proportions (or ancestry coefficients) of each sampled individual.Under the admixture model (Pritchard et al. 2000), an individual i’s ancestry (or admixture proportions) can be characterised by a vector qi = {qi1, qi2, …, qiK}, where qik is the proportion of its genome coming from (contributed by) source population k. Equivalently, qik can also be taken as the probability that an allele sampled at random from individual i comes from source population k. Obviously, we have qik ≥ 0 and (mathop {sum}nolimits_{k = 1}^K {q_{ik} equiv 1}). The overall admixture extent of individual i can be measured by (M_i = 1 – mathop {sum}nolimits_{k = 1}^K {q_{ik}^2}), the probability that the two alleles at a randomly drawn locus come from different source populations. Individual i is purebred and admixed when Mi = 0 and Mi  > 0, respectively. An F1 and F2 hybrid individual i is expected to have Mi = 0.5 and Mi = 0.625, respectively.The task of an admixture analysis is to infer qi for each individual i, denoted by Q = {q1, q2, …, qN}. The log-likelihood function is$${{{mathcal{L}}}}left( {{{{mathbf{Q}}}},{{{mathbf{P}}}}left| {{{mathbf{X}}}} right.} right) = mathop {sum}limits_{i = 1}^N {mathop {sum}limits_{l = 1}^L {mathop {sum}limits_{a = 1}^2 {{{{mathrm{Log}}}}left( {mathop {sum}limits_{k = 1}^K {q_{ik}p_{klx_{ila}}} } right)} } }$$
(6)
Note (6) is essentially the same as those proposed in previous studies (e.g., Tang et al. 2005; Alexander et al. 2009). It assumes independence of individuals conditional on the genetic structure defined by Q, and independence of alleles both within and between loci. The former can be violated when the data have genetic structure in addition to the subpopulation structure defined by Q, such as the presence of familial structure (Rodríguez‐Ramilo and Wang 2012) or inbreeding (Gao et al. 2007) within a subpopulation. The assumption of independence among loci is violated for markers in linkage disequilibrium. It, as well as the assumption of independence between paternal and maternal alleles within a locus, is also violated due to admixture (Tang et al. 2005) or inbreeding (Gao et al. 2007). However, (6) is a good approximation and works well in general even when these assumptions are violated, as checked by extensive simulations.If P were known, it would be trivial to estimate Q from X. Unfortunately, usually, the only information we have is genotype data X, from which we must infer K, Q and P jointly. Herein I modify the EM algorithm of Tang et al. (2005) to solve (6) for maximum likelihood estimates of Q and P given K, as detailed in Supplementary Appendix 2.Despite essentially the same likelihood function, my EM algorithm differs from that of Tang et al. (2005) in three aspects. First, I use the clustering results of mixture model as initial values of Q. Even in the worst scenario of many highly admixed individuals included in a sample, the clustering results should still be much closer to the true Q than a random guess, as used in previous likelihood methods (Tang et al. 2005; Alexander et al. 2009). It is possible (and indeed it has been trialled) to use the results of a faster non-model based clustering method, such as K-means method, in place of those of the likelihood-based clustering method with simulated annealing algorithm as described above. However, such non-model based methods are less reliable and less accurate, especially in difficult situations (below). Second, rather than updating Q and P in alternation, I update Q to asymptotic convergence under a given P. I then update P using the converged Q. This two-step iteration process is repeated until the convergence of both Q and P is reached. Third, the allele frequencies for a specific individual i are calculated by excluding the genotypes of the individual, which are then used in the EM procedure for iteratively updating qi.Optimal KThe above-described clustering analysis and admixture analysis are conducted by assuming a given number of source populations, K. Apparently, different genetic structures would be inferred from the same genotype data if different K values are assumed. In some cases, a reasonable K value is roughly known. For example, individuals might be sampled from K known discrete locations (say, lakes), and the purposes of a structure analysis are to confirm that populations from different locations are indeed differentiated and thus distinguishable, to identify migrants between the locations, and to find out the patterns of genetic differentiations (e.g., whether isolation by distance applies or not). In many other cases, however, we may have no idea of the most likely K value. For example, individuals might be sampled from the same breeding or feeding ground and we wish to know how many populations are using the same ground, and to learn the properties of these populations from the individuals sampled and assigned to them. In such a situation of hidden genetic structure, we need first to identify the most likely one or more K values, and then investigate the corresponding structure/admixture.Estimating the most likely K value from genotype data is difficult (Pritchard et al. 2000). Although many methods have been proposed and applied (see review by Wang 2019), they are all ad hoc to some extent and may be inaccurate in difficult situations such as highly unbalanced sampling from different populations and low differentiation (Wang 2019). Herein I propose two ad hoc estimators of K that can be calculated from the clustering analysis presented in this study. They have a satisfactory accuracy as checked by many test datasets, simulated and empirical.The first estimator is based on the second order rate of change of the estimated log-likelihood as a function of K in a clustering analysis, DLK2. This estimator is similar in spirit to the ∆K method of Evanno et al. (2005), but does not use the mean and standard deviation of log-likelihood values among replicate runs (for a given K value) because the standard deviation (the denominator of ∆K) is frequently zero thanks to the convergence of our clustering analysis by the simulated annealing algorithm.The second estimator, denoted by FSTIS, is based on Wright (1984)’s F-statistics. The best K should produce the strongest population structure, with high differentiation (measured by FST) of each inferred cluster and low deviation from Hardy-Weinberg equilibrium (measured by FIS) within each inferred cluster. Details of how to calculate the two estimators are in Supplementary Appendix 3.SimulationsTo evaluate the accuracy, robustness, and computational efficiency of the new methods implemented in PopCluster in comparison with other methods, I simulated and analysed data with different population structures and sampling intensities. The simulation procedure described below is implemented in the software package PopCluster.Simulation 1, small samplesA population becomes difficult to define genetically when few individuals from it are sampled and included in an admixture analysis. However, a small sample of individuals can be common in practice when, for example, archaeological samples (usually few) are used in studying ancient population structure or in studying the relationship between ancient and current populations (e.g., Lazaridis et al. 2014). In a mixed stock analysis (Smouse et al. 1990) or a wildlife forensic analysis of source populations, there might also be few sampled individuals representing a rare population. To investigate the impact of sample sizes on an admixture analysis, I simulated 10 populations in an island model with FST = 0.05. Nk (=2, 3, …, 10 and 20) individuals were sampled from each of the 10 populations, or 1 individual was sampled from each of the first five populations and 2 individuals were sampled from each of the last five populations (the case Nk = 1.5, Table 1). Other simulation parameters are summarized in Table 1.Table 1 Simulation parameters.Full size tableSimulation 2, many populationsAdmixture becomes increasingly difficult to infer with an increasing K, the number of assumed populations, because the dimensions of both Q and P increase linearly with K. This contrasts with the number of individuals, N, and the number of loci, L, which determines the dimensions of Q and P only, respectively. Therefore, the scale of an admixture analysis, in terms of the number of parameters to be estimated, is predominantly determined by K rather than N or L. I simulated data with a widely variable number of populations (K = [6, 100]) to see if the structure can be accurately reconstructed by using relatively highly informative markers (parameters in Table 1), especially when K is large which is rarely considered in previous simulation studies.Simulation 3, spatial admixture modelThe spatial admixture model resembles isolation by distance where population structure changes gradually as a function of geographic location. Under this model, populations are not discrete as assumed by admixture models and have no recognisable boundaries, posing challenges to an admixture analysis. To simulate the spatially gradual changes in genetic structure, I assume source populations 1, 2, …, K are equally spaced in that order along a line (say, a river in reality). Sampled individuals 1, 2, …, N are also equally spaced in that order on the same line. The admixture proportions of individual i, qi = {qi1, qi2, …, qiK}, being the proportional genetic contributions to i from source populations k, are a function of the individual’s proximity to these K source populations. Formally, we have$$q_{ik} = frac{{q_{ik}^ ast }}{{mathop {sum}nolimits_{k = 1}^K {q_{ik}^ ast } }}$$
(7)
where$$q_{ik}^ ast = left[ {1 – left( {frac{{i – 1}}{{N – 1}} – frac{{k – 1}}{{K – 1}}} right)^2} right]^S$$and parameter S is used to regulate the admixture extent of the N sampled individuals. Under this spatial admixture model, an individual i’s admixture (qi) is determined by its location, or the distances from the K source populations. The 1st and the last sampled individuals (i = 1, N) always have the least admixture, measured by (M_i = 1 – mathop {sum}nolimits_{k = 1}^K {q_{ik}^2}). q11 (=qNK) is always the largest among the qik values for i = 1, 2, …, N and k = 1, 2, …, K. Given a desired value of q11 and K, the scaler parameter S can be solved from the above equations. Given K, N and S, qi of an individual i can then be calculated from the above equations. In this study, I simulated and analysed samples generated with parameters K = 5, N = 500, L = 10000 SNPs, and q11 varying between 0.5 and 1.0 (Table 1).Simulation 4, low differentiationPopulation structure analysis becomes increasingly difficult with a decreasing differentiation, usually measured by FST, among subpopulations. Fortunately, with genomic data of many SNPs, it is still possible to detect weak and subtle population structures (Patterson et al. 2006) as demonstrated in human fine-structure analysis (e.g., Leslie et al. 2015). I simulated data with varying weak population structures (low FST, Table 1) and otherwise ideal populational (only 3 equally differentiated subpopulations) and sampling conditions (i.e., a large sample of individuals per subpopulation, and many SNPs). The number of SNPs used in analyse was L = 1000/FST such that in principle the population structures should be inferred with roughly equal power and accuracy. Because L is large for low FST, STRUCTURE analysis was abandoned due to computational difficulties.Simulation 5, unbalanced samplingSamples of individuals from different source populations are rarely identical in size in practice. Frequently, different source populations are represented by different numbers of individuals in a sample. The impact of unbalanced sampling and how to mitigate it in applying STRUCTURE have been investigated (e.g., Puechmaille 2016; Wang 2017). Similar problems exist for other admixture or clustering analysis methods but have not been studied yet. The same population structure and unbalanced sampling schemes (see parameters in Table 1) used in Wang (2017) were used to simulate data, which were then analysed by various methods to understand their robustness to unbalanced sampling.Simulation 6, computational efficiencySamples from a variable number of populations (Table 1) were analysed by the four programs on a linux cluster to compare their computational efficiencies. Each program uses a single core (no parallelisation) of a processor (Intel Xeon Gold 6248 2.5 GHz) for a maximal allowed time of 48 or 72 (when K = 1024 only) hours. Default parameter settings are used for all four programs. For STRUCTURE, both burn-in and run lengths were set to 104, although much higher burn-in is required for convergence when K is large (say K  > 20). The running time for STRUCTURE is thus conservative, especially when K is not small.Further simulations were conducted to investigate the effects of high admixture and the presence of familial relationships and inbreeding on the relative performance of different admixture analysis methods, as detailed in Supplementary Appendix 4.In all simulations except for the spatial admixture model, I assumed a population with K discrete subpopulations in Wright’s (1931) island model in equilibrium among mutation, drift and migration. For a locus l (=1, 2, …, L) with Jl alleles, allele frequencies of the ancestral population, p0l = {p0l1, p0l2, …, (p_{0lJ_l})}, were drawn from a uniform Dirichlet distribution, ({{{mathcal{D}}}}left( {lambda _1,lambda _2, ldots ,lambda _{J_l}} right)) where λj = 1 for j = 1, 2, …, Jl. Given p0l, allele frequencies of subpopulation k (=1, 2, …, K), pkl = {pkl1, pkl2, …, (p_{klJ_l})}, were drawn from a uniform Dirichlet distribution, ({{{mathcal{D}}}}left( {lambda _1,lambda _2, ldots ,lambda _{J_l}} right)), where (lambda _j = ( {frac{1}{{F_{ST}}} – 1} )p_{0lj}) for j = 1, 2, …, Jl (Nicholson et al. 2002; Falush et al. 2003). Given pkl and the admixture proportion qi of individual i, two alleles at locus l were drawn independently to form the individual’s genotype. The multilocus genotype of an individual was obtained by combining single locus genotypes sampled independently, assuming linkage equilibrium. Nk individuals were drawn at random from population k (= 1, 2, …, K), which were then pooled and subjected to a structure analysis.For the spatial population and sampling model, allele frequencies at a locus l, p0l and pkl, are generated as before, assuming FST = 0.05 among K = 5 subpopulations. A number of N = 500 individuals, equally spaced on the line between source populations 1 and 5, are sampled. The admixture proportion of individual i, qi, is determined by its location, calculated by Eq. (7). Given pkl and qi, the multilocus genotype of individual i is simulated as described above.For each parameter combination, 100 replicate datasets were simulated, analysed and assessed for estimation accuracy. Each dataset was analysed for admixture by different methods (see below for details) with an assumed K as used in simulations. I did not consider estimating the optimal K by analysing a simulated dataset in a range of possible K values. This is because, like previous studies (e.g., Pritchard et al. 2000; Alexander et al. 2009), I am more concerned with admixture inference under a given K, which is important of itself and forms the basis for inferring the optimal K as well. This is also because it is almost impossible computationally to estimate the optimal K for so many replicate datasets and so many parameter combinations in a large-scale simulation study like the present one, even when using large computer clusters. The optimal K was estimated for several empirical datasets (below).Measurement of accuracyInference accuracy could be assessed by comparing, for each individual i, the agreement between simulated ancestry coefficients, qi, and estimated ancestry coefficients, (widehat {{{mathbf{q}}}}_i), obtained by an admixture analysis assuming the true/simulated subpopulation number K. Because the reconstructed populations are labelled arbitrarily (Pritchard et al. 2000), no meaningful results can be gained by comparing qi and (widehat {{{mathbf{q}}}}_i) directly, however. It is possible to relabel the reconstructed populations and find the labelling scheme that has the maximum agreement between qi and (widehat {{{mathbf{q}}}}_i) as the measurement of accuracy. However, there are K! possible labelling schemes, making the approach difficult to calculate when K is large (say, K > 50).The labelling becomes irrelevant when pairs of individuals are considered for the co-assignment probabilities (or coancestry) (Dawson and Belkhir 2001). I calculate and use the average difference between simulated and estimated coancestry for pairs of sampled individuals to measure the average assignment error, AAE (Wang 2017),$$AAE = left( {frac{1}{{Nleft( {N – 1} right)/2}}mathop {sum}limits_{i = 1}^N {mathop {sum}limits_{j = 1 + 1}^N {left( {mathop {sum}limits_{k = 1}^K {q_{ik}q_{jk}} – mathop {sum}limits_{k = 1}^K {widehat q_{ik}widehat q_{jk}} } right)^2} } } right)^{1/2}.$$
(8)
The minimum value of AAE is 0, when ancestry (admixture) is inferred perfectly. The maximum value is 1, when there are no admixed individuals in the sample, individuals from the same source population are always assigned to different populations and individuals from different source populations are always assigned to the same population. It is worth noting that the minimum AAE value of 0 is always possible for any population structure. However, the maximum value varies and can be much smaller than 1, depending on the actual underlying population structure. With an increasing K value or increasing admixture (i.e., qik→1/K for any individual i), the maximum value of AAE tends to decrease. For this reason, AAE cannot be compared fairly between different genetic structures (e.g., different K values, different actual Q for a given K, or different sizes of subsamples from the source populations) for measuring the relative inference qualities. However, it can always be used to compare the accuracy of different inference methods for a given simulated genetic structure and a given sample.Analysis of real datasetsAn ant datasetIt was originally used in a study of the mating system of an ant species, Leptothorax acervorum (Hammond et al. 2001). Ten sampled colonies, A, B, C, D, E, F, G, H, I, and J, contribute respectively 9, 7, 47, 45, 45, 45, 45, 45, 44, and 45 diploid workers to a sample of 377 individuals. For this species, we know that each colony is headed by a single diploid queen mated with a single haploid male. Therefore, workers from the same colony are full-sibs and workers from different colonies are non-sibs. Each sampled worker was genotyped at up to 6 microsatellite loci, which have 3 to 22 alleles per locus observed in the 377 individuals. This dataset was analysed to reconstruct the genetic structure of the sample, which actually is the family structure. ADMIXTURE and sNMF cannot handle multiallelic marker data and therefore only STRUCTURE and PopCluster are used for analysing this dataset.For STRUCTURE, I used the default parameter settings, except for the burning-in and run lengths which were both set to 105 to reduce the risk of non-convergence. Two analyses were conducted. First, optimal K values were determined using three estimators (Wang 2019) calculated from STRUCTURE outputs, and using the DLK2 estimator of PopCluster. For this K estimation purpose, 20 replicate runs for each possible K value in the range [1, 15] were conducted by both STRUCTURE and PopCluster. Second, assuming K = 10, a number of 100 replicate runs (each with a distinctive seed for the random number generator) were conducted by both STRUCTURE and PopCluster to investigate their convergence.An Arctic charr datasetShikano et al. (2015) sampled 328 Arctic charr individuals from 6 locations in northern Fennoscandia: two lakes (Galggojavri and Gallajavri) and one pond (Leenanlampi) in the Skibotn watercourse drain into the Atlantic Ocean and three lakes (Somasjärvi, Urtas-Riimmajärvi and Kilpisjärvi) in the Tornio-Muoniojoki watercourse drain into the Baltic Sea. Individuals were genotyped at 15 microsatellite loci to study the genetic structure and demography. The data were again analysed by STRUCTURE and PopCluster but not by ADMIXTURE and sNMF because the markers are multiallelic. I conducted two separate analyses of the genotype data. First, I estimated the most likely K value by each program, making 20 replicate runs with each K value in the range [1, 10]. Second, I investigated the convergence of each program by conducting 100 replicate runs of the data at K = 6. STRUCTURE analyses were run with default parameter settings except for both burn-in and run lengths being 105.A human SNP datasetUsing FRAPPE (Tang et al. 2005), Li et al. (2008) studied the world-wide human population structure represented by 938 individuals sampled from 51 populations of the Human Genome Diversity Panel (HGDP). Each individual was genotyped at 650000 common SNP loci. The data were expanded to include genotypes of 1043 individuals at 644258 SNPs, available from http://www.cephb.fr/en/hgdp_panel.php#basedonnees. In this study, the expanded data were comparatively analysed by PopCluster, ADMIXTURE, and sNMF, assuming K = 7 clusters (regions) as in the original study (Li et al. 2008). STRUCTURE was too slow to analyse this big dataset and thus it was abandoned.The human 1000 genomes phase I datasetThe dataset (Abecasis et al. 2012), available from https://www.internationalgenome.org/data/, has 1092 human individuals sampled from 14 populations across all continents, with each individual having 38 million SNP genotypes. After removing monomorphic loci (note, no pruning was applied regarding missing data, minor allele frequency and linkage disequilibrium, in contrast to other studies), genotypes at a number of L = 38035992 SNPs were analysed by PopCluster and sNMF, assuming K = 9 clusters (regions). Both STRUCTURE and ADMIXTURE were too slow to analyse this huge dataset and thus were abandoned. No attempts are made to find the optimal K for this dataset as done for the ant and Arctic charr datasets, because too much computational time is required for PopCluster or sNMF to analyse the data with a number of replicate runs at each of a number of K values even when using a large cluster, and there might be multiple K values that explain the data equally well (at different spatial and time scales). For a better understanding of the world-wide human population genetic structure, the data should be analysed at least with one replicate under each of a number of possible K values, say K = [5, 12], to reveal and compare the genetic structure. This study analysed the data at a single K = 9 for the purpose of demonstrating the capacity of different methods, and comparing the admixture estimates of PopCluster and sNMF at this particular value of K. Because of the incompleteness of the analysis, the biological interpretations of the results should be taken with caution.Comparative analyses by different softwareI compared the accuracy and computational time of STRUCTURE (Pritchard et al. 2000; Falush et al. 2003), ADMIXTURE (Alexander et al. 2009), sNMF (Frichot et al. 2014) and PopCluster in analysing both simulated and empirical datasets described above. Quite a few other model-based methods implemented in various software exist. I choose STRUCTURE and ADMIXTURE because they are the most popular model-based admixture analysis methods used for small and large datasets, respectively. I also choose sNMF because it is a very fast model-based method that works for huge datasets for which other methods, such as ADMIXTURE, fail to run or take unrealistically too much time to run.STRUCTURE can handle both diallelic (such as SNPs) and multiallelic (such as microsatellites) markers, but runs too slowly to analyse large datasets with many markers, many individuals, or many populations. It was therefore used to analyse all simulated and empirical datasets with no more than 10000 loci. The default parameter setting was used for most datasets, with a burn-in length of 104 and a run length of 104 iterations. For better convergence, the burn-in and run lengths were increased to 105 iterations for analyses involving a large number of simulated populations (say, when K ≥ 10) or for analyses of empirical datasets. For unbalanced sampling, the alternative ancestry model instead of the default model was used by setting POPALPHAS = 1.Both ADMIXTURE and sNMF were developed specifically for diallelic markers and could not analyse multiallelic marker data. In this study, they were used to analyse SNP data only. For the human 1000 genome phase I data, however, ADMIXTURE could not complete the analysis within a realistic period of time (72 h, the maximum allowed in the linux cluster used for the analysis) even when the maximal number of parallel threads were used. Therefore, only sNMF and PopCluster were used to analyse this dataset.To understand the relative computational efficiency and how much speedup can be gained by parallelisation, ADMIXTURE, sNMF and PopCluster were used to analyse the HGDP dataset and the 1000 genome dataset, by using a variable number of parallel threads on a linux cluster with many nodes, each having 32 cores. The maximum wall clock time allowed for a job on the cluster is 48 h. More

Ollerton, J. Pollinators & Pollination: Nature and Society (Pelagic Publishing, 2021).
Google Scholar 
Delaplane, K. S. Crop pollination By Bees: Evolution, Ecology, Conservation, and Management (CABI, 2021).
Google Scholar 
Aizen, M. A., Garibaldi, L. A., Cunningham, S. A. & Klein, A. M. Long-term global trends in crop yield and production reveal no current pollination shortage but increasing pollinator dependency. Curr. Biol. 18, 1572–1575 (2008).CAS 
PubMed 

Google Scholar 
Koh, I. et al. Modeling the status, trends, and impacts of wild bee abundance in the United States. Proc. Natl. Acad. Sci. 113, 140–145 (2016).ADS 
CAS 
PubMed 

Google Scholar 
Jordan, A., Patch, H. M., Grozinger, C. M. & Khanna, V. Economic dependence and vulnerability of United States agricultural sector on insect-mediated pollination service. Environ. Sci. Technol. 55, 2243–2253 (2021).ADS 
CAS 
PubMed 

Google Scholar 
Biddinger, D. J. & Rajotte, E. G. Integrated pest and pollinator management: Adding a new dimension to an accepted paradigm. Curr. Opin. Insect Sci. 10, 204–209 (2015).PubMed 

Google Scholar 
Egan, P. A., Dicks, L. V., Hokkanen, H. M. T. & Stenberg, J. A. Delivering integrated pest and pollinator management (IPPM). Trends Plant Sci. 25, 577–589 (2020).CAS 
PubMed 

Google Scholar 
Flöhr, A., Stenberg, J. A. & Egan, P. A. The Joint Economic Impact Level (jEIL): A Decision Metric for Integrated Pest and Pollinator Management. In Integrative Biological Control 17–38 (Springer, 2020).
Google Scholar 
Krupke, C. H., Hunt, G. J., Eitzer, B. D., Andino, G. & Given, K. Multiple routes of pesticide exposure for honey bees living near agricultural fields. PLoS ONE 7, e29268 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Traynor, K. S. et al. In-hive pesticide exposome: Assessing risks to migratory honey bees from in-hive pesticide contamination in the Eastern United States. Sci. Rep. 6, 1–16 (2016).
Google Scholar 
Mullin, C. A. et al. High levels of miticides and agrochemicals in North American apiaries: Implications for honey bee health. PLoS ONE 5, e9754 (2010).ADS 
PubMed 
PubMed Central 

Google Scholar 
Ravoet, J., Reybroeck, W. & de Graaf, D. C. Pesticides for apicultural and/or agricultural application found in belgian honey bee wax combs. Bull. Environ. Contam. Toxicol. 94, 543–548 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
Liang, H. C., Bilon, N. & Hay, M. T. Analytical methods for pesticide residues. Water Environ. Res. 86, 2132–2155 (2014).CAS 

Google Scholar 
Calatayud-Vernich, P., Calatayud, F., Simó, E. & Picó, Y. Efficiency of QuEChERS approach for determining 52 pesticide residues in honey and honey bees. MethodsX 3, 452–458 (2016).PubMed 
PubMed Central 

Google Scholar 
Fernandez, M., Pico, Y. & Manes, J. Analytical methods for pesticide residue determination in bee products. J. Food Prot. 65, 1502–1511 (2002).CAS 
PubMed 

Google Scholar 
Vázquez, P. P., Lozano, A., Uclés, S., Ramos, M. M. G. & Fernández-Alba, A. R. A sensitive and efficient method for routine pesticide multiresidue analysis in bee pollen samples using gas and liquid chromatography coupled to tandem mass spectrometry. J. Chromatogr. A 1426, 161–173 (2015).PubMed 

Google Scholar 
Stoner, K. A., Cowles, R. S., Nurse, A. & Eitzer, B. D. Tracking pesticide residues to a plant genus using palynology in pollen trapped from honey bees (Hymenoptera: Apidae) at ornamental plant nurseries. Environ. Entomol. 48, 351–362 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Colwell, M. J., Williams, G. R., Evans, R. C. & Shutler, D. Honey bee-collected pollen in agro-ecosystems reveals diet diversity, diet quality, and pesticide exposure. Ecol. Evol. 7, 7243–7253 (2017).PubMed 
PubMed Central 

Google Scholar 
Favaro, R. et al. Botanical origin of pesticide residues in pollen loads collected by honeybees during and after apple bloom. Front. Physiol. 10, 1069 (2019).PubMed 
PubMed Central 

Google Scholar 
Tosi, S., Costa, C., Vesco, U., Quaglia, G. & Guido, G. A 3-year survey of Italian honey bee-collected pollen reveals widespread contamination by agricultural pesticides. Sci. Total Environ. 615, 208–218 (2018).ADS 
CAS 
PubMed 

Google Scholar 
Chaimanee, V., Chantawannakul, P., Khongphinitbunjong, K., Kamyo, T. & Pettis, J. S. Comparative pesticide exposure to Apis mellifera via honey bee-collected pollen in agricultural and non-agricultural areas of Northern Thailand. J. Apic. Res. 58, 720–729 (2019).
Google Scholar 
Friedle, C., Wallner, K., Rosenkranz, P., Martens, D. & Vetter, W. Pesticide residues in daily bee pollen samples (April–July) from an intensive agricultural region in Southern Germany. Environ. Sci. Pollut. Res. https://doi.org/10.1007/s11356-020-12318-2 (2021).Article 

Google Scholar 
Urbanowicz, C. et al. Low maize pollen collection and low pesticide risk to honey bees in heterogeneous agricultural landscapes. Apidologie 50, 379–390 (2019).
Google Scholar 
Stoner, K. A. & Eitzer, B. D. Using a hazard quotient to evaluate pesticide residues detected in pollen trapped from honey bees (Apis mellifera) in Connecticut. PLoS ONE 8, e77550 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
McArt, S. H., Fersch, A. A., Milano, N. J., Truitt, L. L. & Böröczky, K. High pesticide risk to honey bees despite low focal crop pollen collection during pollination of a mass blooming crop. Sci. Rep. 7, 1–10 (2017).
Google Scholar 
Calatayud-Vernich, P., Calatayud, F., Simó, E., Pascual Aguilar, J. A. & Picó, Y. A two-year monitoring of pesticide hazard in-hive: High honey bee mortality rates during insecticide poisoning episodes in apiaries located near agricultural settings. Chemosphere 232, 471–480 (2019).ADS 
CAS 
PubMed 

Google Scholar 
García-Valcárcel, A. I., Martínez-Ferrer, M. T., Campos-Rivela, J. M. & Hernando Guil, M. D. Analysis of pesticide residues in honeybee (Apis mellifera L.) and in corbicular pollen: Exposure in citrus orchard with an integrated pest management system. Talanta 204, 153–162 (2019).PubMed 

Google Scholar 
Fulton, C. A. et al. An assessment of pesticide exposures and land use of honey bees in Virginia. Chemosphere 222, 489–493 (2019).ADS 
CAS 
PubMed 

Google Scholar 
Calatayud-Vernich, P., Calatayud, F., Simó, E. & Picó, Y. Pesticide residues in honey bees, pollen and beeswax: Assessing beehive exposure. Environ. Pollut. 241, 106–114 (2018).CAS 
PubMed 

Google Scholar 
Ostiguy, N. et al. Honey bee exposure to pesticides: A four-year nationwide study. Insects 10, 13 (2019).PubMed Central 

Google Scholar 
Seeley, T. D. The honey bee colony as a superorganism. Am. Sci. 77, 546–553 (1989).ADS 

Google Scholar 
Thompson, H. M. & Maus, C. The relevance of sublethal effects in honey bee testing for pesticide risk assessment. Pest Manag. Sci. 63, 1058–1061 (2007).CAS 
PubMed 

Google Scholar 
Sponsler, D. B. & Johnson, R. M. Mechanistic modeling of pesticide exposure: The missing keystone of honey bee toxicology. Environ. Toxicol. Chem. 36, 871–881 (2017).CAS 
PubMed 

Google Scholar 
Gradish, A. E. et al. Comparison of pesticide exposure in honey bees (Hymenoptera: Apidae) and Bumble Bees (Hymenoptera: Apidae): implications for risk assessments. Environ. Entomol. 48, 12–21 (2019).PubMed 

Google Scholar 
Tosi, S. & Nieh, J. C. Lethal and sublethal synergistic effects of a new systemic pesticide, flupyradifurone (Sivanto®), on honeybees. Proc. R. Soc. B 286, 20190433 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Iverson, A., Hale, C., Richardson, L., Miller, O. & McArt, S. Synergistic effects of three sterol biosynthesis inhibiting fungicides on the toxicity of a pyrethroid and neonicotinoid insecticide to bumble bees. Apidologie 50, 733 (2019).CAS 

Google Scholar 
Siviter, H. et al. Agrochemicals interact synergistically to increase bee mortality. Nature 596, 389–392 (2021).ADS 
CAS 
PubMed 

Google Scholar 
Graham, K. K. et al. Identities, concentrations, and sources of pesticide exposure in pollen collected by managed bees during blueberry pollination. Sci. Rep. 11, 1–13 (2021).
Google Scholar 
EFSA. Guidance on the risk assessment of plant protection products on bees (Apis mellifera, Bombus spp. and solitary bees). EFSA J. https://doi.org/10.2903/j.efsa.2013.3295 (2013).Article 

Google Scholar 
EPA. Guidance for assessing pesticide risks to bees. (2014).USDA APHIS. Wax Sampling Protocol for the National Honey Bee Disease Survey. (2018).European Committee for Standardization. Foods of plant origin – Multimethod for the determination of pesticide residues using GC- and LC-based analysis following acetonitrile extraction/partitioning and clean-up by dispersive SPE – Modular QuEChERS-method. (2018).Couvillon, M. J. et al. Honey bee foraging distance depends on month and forage type. Apidologie 46, 61–70 (2015).
Google Scholar 
Knight, M. E. et al. An interspecific comparison of foraging range and nest density of four bumblebee (Bombus) species. Mol. Ecol. 14, 1811–1820 (2005).CAS 
PubMed 

Google Scholar 
McArt, S. H., Urbanowicz, C., Mccoshum, S., Irwin, R. E. & Adler, L. S. Landscape predictors of pathogen prevalence and range contractions in US bumblebees. Proc. R. Soc. B 284, 20172181 (2017).PubMed 
PubMed Central 

Google Scholar 
USDA NASS. USDA National Agricultural Statistics Service Cropland Data Layer. (2018).R Core Team. R: A Language and Environment for Statistical Computing. (2019).GraphPad Software. GraphPad Prism. (2017).Hartig, F. DHARMa: Residual Diagnostics for Hierarchical (Multi-Level/Mixed) Regression Models. (2019).Barton, K. MuMIn: Multi-Model Inference. (2019).Fox, J. & Weisburg, S. An {R} Companion to Applied Regression. (2011).Lefcheck, J. S. piecewiseSEM: Piecewise structural equation modeling in R for ecology, evolution, and systematics. Methods Ecol. Evol. 7, 573–579 (2015).
Google Scholar 
Hothorn, T., Bretz, F. & Westfall, P. Simultaneous inference in general parametric models. Biom. J. 50, 346–363 (2008).MathSciNet 
PubMed 
MATH 

Google Scholar 
Fox, J. RcmdrMisc: R Commander Miscellaneous Functions. (2020).Bhattacharya, M., Primack, R. B. & Gerwein, J. Are roads and railroads barriers to bumblebee movement in a temperate suburban conservation area?. Biol. Conserv. 109, 37–45 (2003).
Google Scholar 
Fragoso, F. P. & Brunet, J. Patch fidelity of honey bees and bumble bees differs and is affected by spatial configuration. In Entomological Society of America Annual Meeting, Plant-Insect Ecosystems (2021).Javorek, S. K., Mackenzie, K. E. & Vander Kloet, S. P. Comparative (Hymenoptera: Apoidea) on Lowbush Blueberry (Ericaceae: Vaccinium angustifolium). Ann. Entomol. Soc. Am. 95, 345–351 (2002).
Google Scholar 
Sandrock, C. et al. Impact of chronic neonicotinoid exposure on honeybee colony performance and queen supersedure. PLoS ONE 9, e103592 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Wood, S. C., Kozii, I. V., Koziy, R. V., Epp, T. & Simko, E. Comparative chronic toxicity of three neonicotinoids on New Zealand packaged honey bees. PLoS ONE 13, e0190517 (2018).PubMed 
PubMed Central 

Google Scholar 
Degrandi-Hoffman, G., Graham, H., Ahumada, F., Smart, M. & Ziolkowski, N. The economics of honey bee (Hymenoptera: Apidae) management and overwintering strategies for colonies used to pollinate almonds. J. Econ. Entomol. https://doi.org/10.1093/jee/toz213 (2019).Article 
PubMed 

Google Scholar 
Biddinger, D. J. et al. Comparative toxicities and synergism of apple orchard pesticides to Apis mellifera (L.) and Osmia cornifrons (Radoszkowski). PLoS ONE 8, e72587 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Johnson, R. M., Dahlgren, L., Siegfried, B. D. & Ellis, M. D. Acaricide, fungicide and drug Interactions in Honey Bees (Apis mellifera). PLoS ONE 8, e54092 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Jaffe, B. D., Lois, A. N. & Guédot, C. Effect of fungicide on pollen foraging by honeybees (Hymenoptera: Apidae) in cranberry differs by fungicide type. J. Econ. Entomol. 112, 499–503 (2019).CAS 
PubMed 

Google Scholar 
Sanchez-Bayo, F. & Goka, K. Pesticide residues and bees: A risk assessment. PLoS ONE 9, e94482 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Wang, Y., Zhu, Y. C. & Li, W. Interaction patterns and combined toxic effects of acetamiprid in combination with seven pesticides on honey bee (Apis mellifera L.). Ecotoxicol. Environ. Saf. 190, 110100 (2020).CAS 
PubMed 

Google Scholar 
Wang, Y., Zhu, Y. C. & Li, W. Comparative examination on synergistic toxicities of chlorpyrifos, acephate, or tetraconazole mixed with pyrethroid insecticides to honey bees (Apis mellifera L.). Environ. Sci. Pollut. Res. 27, 6971–6980 (2019).
Google Scholar 
Sgolastra, F. et al. Synergistic mortality between a neonicotinoid insecticide and an ergosterol-biosynthesis-inhibiting fungicide in three bee species. Pest Manag. Sci. 73, 1236–1243 (2017).CAS 
PubMed 

Google Scholar 
Azpiazu, C. et al. Chronic oral exposure to field-realistic pesticide combinations via pollen and nectar: effects on feeding and thermal performance in a solitary bee. Sci. Reports 9, 1–11 (2019).CAS 

Google Scholar 
Becher, M. A., Hildenbrandt, H., Hemelrijk, C. K. & Moritz, R. F. A. Brood temperature, task division and colony survival in honeybees: A model. Ecol. Modell. 221, 769–776 (2010).
Google Scholar 
Zhu, W., Schmehl, D. R., Mullin, C. A. & Frazier, J. L. Four common pesticides, their mixtures and a formulation solvent in the hive environment have high oral toxicity to honey bee larvae. PLoS ONE 9, e77547 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Dai, P. et al. Chronic toxicity of clothianidin, imidacloprid, chlorpyrifos, and dimethoate to Apis mellifera L. larvae reared in vitro. Pest Manag. Sci. 75, 29–36 (2019).CAS 
PubMed 

Google Scholar 
Campbell, J. B. et al. The fungicide Pristine® inhibits mitochondrial function in vitro but not flight metabolic rates in honey bees. J. Insect Physiol. 86, 11–16 (2016).CAS 
PubMed 

Google Scholar 
DesJardins, N. S. et al. A common fungicide, Pristine®, impairs olfactory associative learning performance in honey bees (Apis mellifera). Environ. Pollut. 288, 117720 (2021).CAS 
PubMed 

Google Scholar 
Fisher, A. et al. Colony field test reveals dramatically higher toxicity of a widely-used mito-toxic fungicide on honey bees (Apis mellifera). Environ. Pollut. 269, 115964 (2021).CAS 
PubMed 

Google Scholar 
Mahefarisoa, K. L., Simon Delso, N., Zaninotto, V., Colin, M. E. & Bonmatin, J. M. The threat of veterinary medicinal products and biocides on pollinators: A one health perspective. One Heal. 12, 100237 (2021).CAS 

Google Scholar 
Christen, V., Schirrmann, M., Frey, J. E. & Fent, K. Global transcriptomic effects of environmentally relevant concentrations of the neonicotinoids clothianidin, imidacloprid, and thiamethoxam in the brain of honey bees (Apis mellifera). Environ. Sci. Technol. 52, 7534–7544 (2018).ADS 
CAS 
PubMed 

Google Scholar 
Tison, L., Duer, A., Púčiková, V., Greggers, U. & Menzel, R. Detrimental effects of clothianidin on foraging and dance communication in honey bees. PLoS ONE 15, e0241134 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Tison, L., Rößner, A., Gerschewski, S. & Menzel, R. The neonicotinoid clothianidin impairs memory processing in honey bees. Ecotoxicol. Environ. Saf. 180, 139–145 (2019).CAS 
PubMed 

Google Scholar 
Morfin, N., Goodwin, P. H., Correa-Benitez, A. & Guzman-Novoa, E. Sublethal exposure to clothianidin during the larval stage causes long-term impairment of hygienic and foraging behaviours of honey bees. Apidologie 50, 595–605 (2019).CAS 

Google Scholar 
Yao, J., Zhu, Y. C. & Adamczyk, J. Responses of honey bees to lethal and sublethal doses of formulated clothianidin alone and mixtures. J. Econ. Entomol. 111, 1517–1525 (2018).CAS 
PubMed 

Google Scholar 
Bortolotti, L. et al. Effects of sub-lethal imidacloprid doses on the homing rate and foraging activity of honey bees. Bull. Insectology 56, 63–67 (2003).
Google Scholar 
Yang, E. C., Chuang, Y. C., Chen, Y. L. & Chang, L. H. Abnormal foraging behavior induced by sublethal dosage of imidacloprid in the honey bee (Hymenoptera: Apidae). J. Econ. Entomol. 101, 1743–1748 (2008).CAS 
PubMed 

Google Scholar 
Karahan, A., Cakmak, I., Hranitz, J. M., Karaca, I. & Wells, H. Sublethal imidacloprid effects on honey bee flower choices when foraging. Ecotoxicology 24, 2017–2025 (2015).CAS 
PubMed 

Google Scholar 
Dively, G. P., Embrey, M. S., Kamel, A., Hawthorne, D. J. & Pettis, J. S. Assessment of chronic sublethal effects of imidacloprid on honey bee colony health. PLoS ONE 10, 1–25 (2015).
Google Scholar 
Meikle, W. G. et al. Sublethal effects of imidacloprid on honey bee colony growth and activity at three sites in the U.S.. PLoS ONE https://doi.org/10.1371/journal.pone.0168603 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Wu, Y.-Y. et al. Sublethal effects of imidacloprid on targeting muscle and ribosomal protein related genes in the honey bee Apis mellifera L.. Sci. Rep. https://doi.org/10.1038/s41598-017-16245-0 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Wu, M.-C., Chang, Y.-W., Lu, K.-H. & Yang, E.-C. Gene expression changes in honey bees induced by sublethal imidacloprid exposure during the larval stage. Insect Biochem. Mol. Biol. 88, 12–20 (2017).CAS 
PubMed 

Google Scholar 
Peng, Y.-C. & Yang, E.-C. Sublethal dosage of imidacloprid reduces the microglomerular density of honey bee mushroom bodies. Sci. Rep. https://doi.org/10.1038/srep19298 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
de Almeida-Rossi, C., Roat, T. C., Tavares, D. A., Cintra-Socolowski, P. & Malaspina, O. Brain morphophysiology of africanized bee Apis mellifera exposed to sublethal doses of imidacloprid. Arch. Environ. Contam. Toxicol. 65, 234–243 (2013).PubMed 

Google Scholar 
Tosi, S., Burgio, G. & Nieh, J. C. A common neonicotinoid pesticide, thiamethoxam, impairs honey bee flight ability. Sci. Reports 7, 1–8 (2017).
Google Scholar 
Coulon, M. et al. Interactions between thiamethoxam and deformed wing virus can drastically impair flight behavior of honey bees. Front. Microbiol. 0, 766 (2020).
Google Scholar 
Shi, T.-F., Wang, Y.-F., Liu, F., Qi, L. & Yu, L.-S. Sublethal effects of the neonicotinoid insecticide thiamethoxam on the transcriptome of the honey bees (Hymenoptera: Apidae). J. Econ. Entomol. 110, 2283–2289 (2017).CAS 
PubMed 

Google Scholar 
Tesovnik, T. et al. Exposure of honey bee larvae to thiamethoxam and its interaction with Nosema ceranae infection in adult honey bees. Environ. Pollut. 256, 113443 (2020).CAS 
PubMed 

Google Scholar 
Friol, P. S., Catae, A. F., Tavares, D. A., Malaspina, O. & Roat, T. C. Can the exposure of Apis mellifera (Hymenoptera, Apiadae) larvae to a field concentration of thiamethoxam affect newly emerged bees?. Chemosphere 185, 56–66 (2017).ADS 
CAS 
PubMed 

Google Scholar 
Coulon, M. et al. Influence of chronic exposure to thiamethoxam and chronic bee paralysis virus on winter honey bees. PLoS ONE 14, e0220703 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Stoner, K. A. & Eitzer, B. D. Movement of soil-applied imidacloprid and thiamethoxam into nectar and pollen of squash (Cucurbita pepo). PLoS ONE 7, e39114 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Rolke, D., Persigehl, M., Peters, B., Sterk, G. & Blenau, W. Large-scale monitoring of effects of clothianidin-dressed oilseed rape seeds on pollinating insects in northern Germany: residues of clothianidin in pollen, nectar and honey. Ecotoxicology 25, 1691–1701 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Wood, T. J., Kaplan, I., Zhang, Y. & Szendrei, Z. Honeybee dietary neonicotinoid exposure is associated with pollen collection from agricultural weeds. Proc. R. Soc. B 286, 1905 (2019).
Google Scholar 
Urlacher, E. et al. Measurements of chlorpyrifos levels in forager bees and comparison with levels that disrupt honey bee odor-mediated learning under laboratory conditions. J. Chem. Ecol. 42, 127–138 (2016).CAS 
PubMed 

Google Scholar 
Villalba, A., Maggi, M., Ondarza, P. M., Szawarski, N. & Miglioranza, K. S. B. Influence of land use on chlorpyrifos and persistent organic pollutant levels in honey bees, bee bread and honey: Beehive exposure assessment. Sci. Total Environ. 713, 136554 (2020).ADS 
CAS 
PubMed 

Google Scholar 
Traynor, K. S. et al. Pesticides in honey bee colonies: Establishing a baseline for real world exposure over seven years in the USA. Environ. Pollut. https://doi.org/10.1016/j.envpol.2021.116566 (2021).Article 
PubMed 

Google Scholar 
EPA Press Office. EPA Takes Action to Address Risk from Chlorpyrifos and Protect Children’s Health. (2021).Arena, M. & Sgolastra, F. A meta-analysis comparing the sensitivity of bees to pesticides. Ecotoxicology 233(23), 324–334 (2014).
Google Scholar 
Wright, G. A., Nicolson, S. W. & Shafir, S. Nutritional physiology and ecology of honey bees. Annu. Rev. Entomol 63, 327–344. https://doi.org/10.1146/annurev-ento-020117-043423 (2018).CAS 
Article 
PubMed 

Google Scholar 
Michener, C. D. The Social Behavior of the Bees: A Comparative Study (Belknap Press of Harvard University Press, 1974).
Google Scholar 
Stabler, D., Paoli, P. P., Nicolson, S. W. & Wright, G. A. Nutrient balancing of the adult worker bumblebee (Bombus terrestris) depends on the dietary source of essential amino acids. J. Exp. Biol. https://doi.org/10.1242/jeb.114249 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Böhme, F., Bischoff, G., Zebitz, C. P. W., Rosenkranz, P. & Wallner, K. Pesticide residue survey of pollen loads collected by honeybees (Apis mellifera) in daily intervals at three agricultural sites in South Germany. PLoS ONE 13, e0199995 (2018).PubMed 
PubMed Central 

Google Scholar 
Schilder, A. M. C., Hanson, E. J. & Hancock, J. F. An integrated approach to disease control in blueberries in Michigan. Acta Hortic. 715, 481–488 (2006).CAS 

Google Scholar 
Wise, J. C., Jenkins, P. E., Poppen, R. V. & Isaacs, R. Activity of broad-spectrum and reduced-risk insecticides on various life stages of cranberry fruitworm (Lepidoptera: Pyralidae) in Highbush Blueberry. J. Econ. Entomol. 103, 1720–1728 (2010).PubMed 

Google Scholar 
Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2009).MATH 

Google Scholar  More

Study sites, sample collection and preparationAdult female mosquitoes used in this study had previously been collected from three areas: Kerio Valley (Baringo county), Rabai (Kilifi county) and Nguruman (Kajiado county) (Fig. 1), as part of vector-borne disease surveillance project and stored at – 80 °C at the International Centre of Insect Physiology and Ecology (icipe). The mosquitoes were surveyed between August 2019 and May 2020. Nguruman is an agropastoral area located in Kajiado county at the southern end of the Kenyan Rift Valley bordering Tanzania. The area has a semi-arid climate characterized by erratic rains, extreme temperatures, and cyclic and prolonged droughts30. The vegetation is dominated by bushland, grassland and open woodlands along seasonal river valleys. Specific indicator data for malaria is not available for Nguruman except for estimates pertaining to the larger Kajiado county which as of 2019 indicates a malaria incidence rate of 5 per 1000 population31. Collections in Kerio Valley (Baringo county within the Rift Valley) were conducted in Kapluk and Barwesa, both agro-pastoral areas with arid and semi-arid ecology. Malaria is a major vector-borne disease in the areas with report of perennially occurrence in neighboring riverine areas32. Rabai is one of the seven administrative sub-counties of Kilifi county in the coastal region of Kenya where malaria is endemic. The main economic activities in the area include subsistence agriculture, casual labor, crafts and petty trading. The weather patterns at the sites during the sampling period were as follows: Kerio Valley (mean daily temperature: 21.2 °C, mean daily rainfall: 4.1 mm, mean relative humidity: 73.4%); Rabai (mean daily temperature: 26.4 °C, mean daily rainfall: 2.1 mm; mean relative humidity: 78.1%) and Nguruman (mean daily temperature: 22.5 °C, mean daily rainfall: 0.9 mm, mean relative humidity: 61.2%).Mosquito survey and processingHost seeking mosquitoes were trapped using CDC light traps baited with dry ice (carbon dioxide) attractive to several mosquitoes. Traps were set outdoors about 10–15 m away from randomly selected homesteads from 18:00 h to 06:00 h. After collection, the mosquitoes were anesthetized with trimethylamine and temporarily stored in liquid nitrogen before transportation to the Emerging Infectious Disease (EID) laboratory at icipe and later stored at − 80 °C. Anopheline mosquitoes were morphologically identified to species level using published taxonomic keys15,33.DNA extraction and Anopheles species discriminationDNA was extracted from the head/thorax of individual mosquitoes using ISOLATE II Genomic DNA Extraction kit (Bioline, UK) following the manufacturer’s instructions and used for species discrimination and screening for P. falciparum infection and Gste2 mutations (described below).Cryptic sibling species of the Anopheles funestus and Anopheles gambiae complexes were identified using conventional PCR34,35 and/or sequencing. PCR for An. funestus complex in a 15 µl reaction volume comprised 0.5 µM of each primer targeting: Anopheles funestus s.s, Anopheles vaneedeni, Anopheles rivulorum, Anopheles parensis, Anopheles leesoni, Anopheles longipalpis A and Anopheles longipalpis C, 3 µl of 5X HOT FIREPol Blend Master Mix Ready to Load (Solis BioDyne, Estonia) and 2 µl of DNA template. The cycling conditions were initial denaturation at 95 °C for 15 min, and then 30 cycles of denaturation at 95 °C for 30 s, annealing at 46 °C for 30 s and extension at 72 °C for 40 s and final extension at 72 °C for 10 min. Size fragments of each species were scored after separation in 1.5% agarose gel electrophoresis stained with ethidium bromide against a 1 Kb DNA ladder (HyperLadder, Bioline, London, UK).For An. gambiae s.l., PCR in a 10 µl volume consisted of 2 µl of 5X Evagreen HRM Master Mix (Solis BioDyne, Estonia), 1 µl of DNA template and 10 µM concentration of each primer targeting An. gambiae s.s and An. arabiensis. The thermal cycling conditions included initial denaturation for 15 min at 95 °C followed by 40 cycles of denaturation at 95 °C for 20 s, annealing at 61 °C for 15 s and extension at 72 °C for 20 s followed by final extension at 72 °C for 7 min.A subset of An. funestus s.l. samples that failed to amplify using the established protocol, was further amplified and sequenced targeting the internal transcribed spacer 2 (ITS2) region of the ribosomal DNA (rDNA)36. This target has shown utility in discriminating closely related mosquito species including anophelines12 and sequences from diverse species for this marker are well represented in reference databases (e.g. GenBank). PCR volumes for rDNA ITS2 were 15 µl containing 0.5 µM of the forward and reverse primers, 3 µl of 5X HOT FIREPol Blend Master Mix Ready to Load (Solis BioDyne, Estonia) and 2 µl of DNA template. The cycling conditions were initial denaturation at 95 °C for 15 min, followed by 40 cycles of denaturation at 95 °C for 30 s, annealing at 60 °C for 30 s and extension at 72 °C for 45 s and final extension at 72 °C for 7 min. ExoSAP IT rapid cleanup kit (Affymetrix Inc., Santa Clara, CA, USA) was used to clean the PCR product as per the manufacturer’s guideline, and then outsourced for bidirectional Sanger sequencing to Macrogen, South Korea.Detection of malaria parasitesPlasmodium falciparum sporozoites in individual mosquitoes (head/thorax) were detected by analyzing high resolution melting (HRM) profiles generated from real time PCR products of non-coding mitochondrial sequence (ncMS)37. A P. falciparum DNA from National Institute for Biological Standards and Control (NIBSC; London, UK) was used as a reference positive control. PCR was carried out in a 10 µl volume consisting of 2 µl of 5X Evagreen HRM Master Mix (Solis BioDyne, Estonia), 1 µl of DNA template and 10 µM of each primer. PCR cycling conditions were initial denaturation for 15 min at 95 °C followed by 40 cycles of denaturation at 95 °C for 20 s, annealing at 61 °C for 15 s and extension at 72 °C for 20 s followed by final extension at 72 °C for 7 min. A fraction of RT-PCR-HRM positive samples were further analyzed using conventional PCR in a 10 µl volume consisting of 2 µl of 5X HOT FIREPol Blend Master Mix Ready to Load (Solis BioDyne, Estonia), 1 µl of DNA template and 10 µM of each primer. The cycling conditions comprised initial denaturation for 15 min at 95 °C followed by 40 cycles of denaturation at 95 °C for 20 s, annealing at 61 °C for 15 s and extension at 72 °C for 20 s followed by final extension at 72 °C for 7 min. PCR product of samples positive by RT-PCR were purified using ExoSAP- IT (USB Corporation, Cleveland, OH, USA) and outsourced for sequencing to Macrogen, South Korea. All sporozoite-positive mosquitoes were molecularly identified to species by PCR of the ITS2 region as described above.Genotyping for L119F-GSTe2 mutation and sequencingTwo outer and two inner primers in a PCR assay were used to genotype the L119F-GSTe2 mutations that confer resistance of An. funestus mosquitoes to pyrethroids/DDT19 as described previously28. Thus, only An. funestus s.l. was screened using this assay. Briefly, PCR in a 15 µl reaction volume consisted of 10 µM of each primer, 3 µl of 5X HOT FIREPol Blend Master Mix Ready to Load (Solis BioDyne, Estonia), and 2 µl of DNA template. The cycling conditions were initial denaturation at 95 °C for 15 min, followed by 40 cycles of denaturation at 95 °C for 30 s, annealing at 59 °C for 30 s and extension at 72 °C for 40 s and final extension at 72 °C for 7 min. Amplicons were resolved in a 1.5% agarose gel stained with ethidium bromide (Sigma-Aldrich, GmbH, Germany) against a 1 Kb DNA ladder (HyperLadder, Bioline, London, UK). The amplicons were scored as either homozygous susceptible (SS) at 312 bp, homozygous resistant (RR) at 523 bp or heterozygous (RS) when both bands were visualized.Representative GSTe2 allele positive samples were sequenced for the GSTe2 gene using the Gste2F and Gste2R primers as described previously38. PCR comprised a reaction volume of 15 µl in MyTaq DNA Polymerase Kit (Bioline, London, UK) containing 10 µM of each primer, 5X My Taq reaction buffer, 2 µl of My taq DNA polymerase and 1 µl of DNA template. PCR conditions were: initial denaturation of 5 min at 95 °C, followed by 30 cycles of 94 °C for 30 s, 58 °C for 30 s and 72 °C for 1 min, with a final extension at 72 °C for 10 min. Cleaning and sequencing of amplicons were performed as described above.Sequence and polymorphism analysisSequences (mosquito, P. falciparum, GSTe2) were viewed and cleaned in Geneious Prime39 and queried in GenBank using Basic Local Alignment Search Tool (BLastn). Parasite sequences were assigned as P. falciparum after  > 98% percentage identity. MAFFT in Geneious Prime39 was used to perform multiple sequence alignments with default parameters. Maximum likelihood (ML) trees were inferred for mosquito ITS2 sequences using the best fit model of sequence evolution with nodal support for different groupings evaluated through 1000 bootstrap replications. GSTe2 gene polymorphism analysis was performed in Geneious Prime and ML tree reconstructed from MAFFT alignment using PhyML v. 2.2.4. Haplotype distribution network was constructed using Templeton-Crandall Sing (TCS) program v. 1.2140.Statistical analysisRelative abundance was used to estimate the outdoor composition of the anopheline mosquitoes. Daily counts of female mosquito/trap/night for An. funestus s.l. and An. gambiae s.l. were compared for each area using generalized linear models (GLM) with negative binomial error structure based on best-fit model residuals. The mean catches/trap/night was computed for each of the species complexes. The P. falciparum sporozoite infection rates (Pfsp) were expressed as the number of positive specimens of the total number of specimens examined. The distribution of L119F-GSTe2 mutations was assessed by determining allelic frequencies in different species. Infection status among the resistant mosquitoes was compared using the Fisher’s Exact Test. Data were analyzed using R v 4.1.0 software at 95% confidence limit.Ethical considerationsEthical review and approval of the study was granted by the Scientific and Ethical Review Unit (SERU) of the Kenya Medical Research Institute (KEMRI) (Protocol No. SSC 2787). Prior to data collection, the purpose of the study, procedures and associated benefits/risks were provided to the local leadership at county and community levels. Additionally, informed verbal consent to trap mosquitoes around homesteads was obtained from household heads. More

Argonema gen. nov. Skoupý et Dvořák.Type species: Argonema galeatum.Morphology: Filamentous cyanobacterium, colonies macroscopic, growing in round bulbs and tufts. The filaments are dark green to blue-green, grey-green or brown-green in color. Cells are wider than they are long. Filaments sheathed, sheaths are colorless to light brown, distinct, and variable in length. The filament can protrude from the sheath or the sheath can exceed filament. Trichomes are cylindrical, not attenuated to slightly attenuated towards the end, slightly or not constricted at cell walls. The apical cell can be concave, dark brown, purple-brown to almost black. Cell content often granulated. Necridic cells present, reproduction by hormogonia. The morphological description was based on both culture and fresh material.Etymology: The genus epithet (Argonema) is derived from greek Argo – slow, latent (αργός) and nema – thread (νήμα).A. galeatum sp. nov. Skoupý et Dvořák.Morphology: The cells of A. galeatum are 6.5–9.1 µm (mean 7.81 µm) wide and 1.1–2.5 µm (mean 1.83 µm) long (Figs. 1–5). Filaments are straight, blue-green to gray-green in color. The sheaths are colorless to light brown, distinct, and variable in length. The filament can protrude from the sheath or the sheath can exceed filament. No true branching was observed. Trichomes are cylindrical, not attenuated or slightly attenuated towards the end, slightly or not constricted at cell walls. Some filaments have a concave apical cell that is dark brown, purple-brown to almost black (Fig. 11b). Cell content often granulated. Reproduction by necridic cells and subsequent breaking of the filaments into hormogonia (Fig. 11a,c). The morphological description was based on both culture and fresh material.Figures 1-8Microphotographs of Argonema galeatum (Figs 1–5) and Argonema antarcticum (Figs. 6–8) Trichomes of A. galeatum appear more straight (Fig 2), while trichomes of A. antarcticum form waves (Fig 6) and loops (Fig 7). Scale = 10 µm, wide arrow = necridic cells, arrowhead = granules, asterisk = colored apical cell, circle = empty sheath.Full size imageFigures 9 and 10Histograms of cell dimensions constructed using PAST software. Fig. 9 – Histogram of cell width frequencies in A. galeatum (blue) and A. antarcticum (red). Fig. 10 – Histogram of cell length frequencies in A. galeatum (blue) and A. antarcticum (red).Full size imageHolotype: 38,057, Herbarium of the Department of Botany (OL), Palacký University Olomouc, Czech Republic.Reference strain: Argonema galeatum A003/A1.Type locality: James Ross Island, Western Antarctica, 63.80589S, 57.92147 W.Habitat: Well-developed soil crust.Etymology: Species epithet A. galeatum was derived from latin galea – helmet.A. antarcticum sp. nov. Skoupý et Dvořák.Morphology: The cells are 7.6–9.2 µm (mean 8.52 µm) wide and 1.2–2.8 µm (mean 1.72 µm) long (Figs. 5–8). Filaments are wavy, gray-green to brown-green in color. The sheaths are colorless to light brown, distinct, and variable in length. The filament can protrude from the sheath or the sheath can exceed filament. No true branching was observed. Trichomes are cylindrical, not attenuated or slightly attenuated towards the end with a concave apical cell, slightly or not constricted at cell walls (Fig. 11d). Necridic cells present (Fig. 11e), reproduction by hormogonia. The morphological description was based on both culture and fresh material.Holotype: 38,058, Herbarium of the Department of Botany (OL), Palacký University, Olomouc, Czech Republic.Reference strain: Argonema antarcticum A004/B2.Type locality: James Ross Island, Western Antarctica, 63.89762S, 57.79743 W.Habitat: Well-developed soil crust.Etymology: Species epithet A. antarcticum was derived from the original sampling site.Morphological variabilityWe used light microscopy to assess the morphology of Argonema from soil crust samples and cultured strains. Argonema is morphologically similar to other Oscillatoriales, such as Lyngbya, Phormidium, and Oscillatoria. In culture, the morphology of A. galeatum and A. antarcticum differed slightly. Filaments of A. antarcticum are wider than cells of A. galeatum, averaging at 8.52 µm (A. galeatum – 7.81 µm). The average cell width/length ratio is 4.54 for A.galeatum and 4.89 for A. antarcticum. The cell width was significantly different between the two species (Nested ANOVA, p  More