Philippa Kaur

Influence of the PLHP on water depth distributionThe hydrodynamic process of Poyang Lake with and without the PLHP is simulated by M1 and M2, respectively. By comparing and analyzing the simulation results, we obtain the changes of the water depth, i.e., the water depth in M2 minus the water depth in M1, in Poyang Lake. Figure 5 shows the mean monthly water depth differences during September and October in three typical years. The water depth in Poyang Lake has increased obviously in most cases after the operation of the project. Combining the water level processes in Fig. 4 and the water depth differences in Fig. 5, the area of the changes in the water depth is seen to be mainly controlled by the Higher Water Levels (HWL) between M1 and M2. While the magnitude of the changes is mainly controlled by the Differences in Water Levels (DWL) between M1 and M2.Figure 5Water depth variation during September and October in the typical years. The figure was generated by Tecplot2020 (https://www.tecplot.com/).Full size imageAs shown in the water level variations processes in Fig. 4a, in the low-water-level year (2006), where the natural inflow is relatively small, the water level in M2 is much higher than that in M1 and reaches the peak value around 10th of October. The DWL also reaches maximum in this time and then gradually decreases. As a consequence, the water depth increases the most in October 2006. As shown in Fig. 5d, in most areas of Poyang Lake, the Xiuhe River and the Ganjiang River, the water depth is increased by more than 1 m. Especially in the main channel from the PLHP to Tangyin and Wucheng, the increase can exceed 4 m, as shown in the red parts in Fig. 5d. While over the surrounding flooded land, the increase ranges from 1 to 3 m. The increase in the water depth during September was essentially the same as that in October, but the area and magnitude of the changes are slightly reduced. The maximum increase is about 3.1 m, which is mainly concentrated in the main channel on the north of Songmen Mountain.In the medium-water-level year (2018), the water level in M1 is constantly lower than that in M2 during September and October, which can be seen in Fig. 4b. Both the HWL and DWL first increase and then decrease, reaching their maximum values in early September and late September, respectively. Furthermore, the mean monthly HWL and DWL in September are slightly greater than those in October. As a consequence, the area and magnitude of the changes in Fig. 5b are larger than those in Fig. 5e. In September 2018, almost the entire main lake region becomes influenced by the PLHP and the increase in water depth ranges from 0 to 3.4 m. While in October 2018, the increase in water depth ranges from 0 to 3.2 m, and the significant increase is mainly concentrated in the main channel on the north of Songmen Mountain.In the high-water-level year (2010), as shown in Fig. 4c, the HWL rises and the DWL decreases as compared to those in other years. During the first 35 days, the hydrodynamic conditions in M2 are exactly the same as that in M1, so the PLHP has no effect on the water depth in the lake region in September (Fig. 5c). However, after the 6th of October, the water level in M2 begins to be higher than that in M1 under the regulation of the PLHP. Although the DWL is small during this time, the HWL is relatively high and thus results in large areas of water depth increase, as shown in Fig. 5f. The maximum increase in water depth is approximately 1.1 m.In these two months, the northeastern parts of the two national nature reserves are observably influenced., The area of the changes is greatest in September 2018, there has been a marked increase in water depth in most areas of the reserves, with the maximum increase reaching about 2.8 m. While in other periods during September and October, with the exception of September 2010, the water depth is increased significantly in about 1/3 ~ 1/2 area of the reserves, with the maximum increase reaching about 2.6 m. In the meantime, the water depth in the southwestern part of these two reserves remains essentially the same as before the operation of the PLHP, with an increase of less than 0.25 m as shown in the grey parts in Fig. 5.Influence of the PLHP on habitat suitability of Vallisneria natansAccording to the relationship between the habitat suitability of Vallisneria natans and the water depth as mentioned in Fig. 3, the water depth results can be translated into the habitat suitability of Vallisneria natans, as shown in Fig. 6. The first and third rows are the distributions of the habitat suitability during September and October, respectively, in the three typical years before the operation of the PLHP, while the second and fourth rows are distributions of the habitat suitability after the operation of the PLHP. The grey parts in Fig. 6 indicate that the habitat suitability is 0, implying that the area is dry or the water depth is greater than 4 m.Figure 6Distributions of the habitat suitability of Vallisneria natans in typical years without (the first and third rows) and with (the second and fourth rows) the PLHP. The figure was generated by Tecplot2020 (https://www.tecplot.com/).Full size imageAs shown in Fig. 6, the suitable area for the growth of Vallisneria natans is mainly concentrated over the flooded land. As the water depth in the main channel is generally more than 4 m, so the habitat suitability is usually 0 there.Before the operation of the PLHP (M1), the water level in 2010 is higher than that those in 2006 and 2018. Therefore, there are more areas covered by water in 2010, and the habitat suitability in Fig. 6c is greater than those in Fig. 6a,b. Similarly, the habitat suitability in Fig. 6i is greater than those in Fig. 6g,h. Because the lake bed is lower in northeastern part than that in the southwestern part, the water depth generally decreases from northeast to southwest. During September and October in 2006 and 2018, large areas of bed in the northeastern part are covered by water, and the water depth is less than 4 m, which is suitable for the growth of Vallisneria natans. While large areas in the southwestern part of the lake are dry, and thus the habitat suitability is 0, as shown in the grey parts in the southwestern part of the main lake region in Fig. 6a,b,g,h. However, because the water level is relatively high during September and October in 2010, there are almost no dry areas in the lake region. As a result, with the exception of the main channel, where the water depth is greater than 4 m, most of the areas in the lake region are suitable for the growth of Vallisneria natans. The most suitable areas for Vallisneria natans vary between these two months, as the water depth in September 2010 is greater than that in October 2010. As shown in Fig. 6c, the red parts are mainly concentrated in the southwestern part of the lake, because there are large areas of flooded region with water depths ranging from 1 to 2 m, which is ideal for the growth of Vallisneria natans. On the contrary, the water depth in the northeastern flood land is usually between 2-4 m. In Fig. 6i, however, the red parts are mainly concentrated in the northeastern part of the lake, because the water depth there is usually between 1 and 2 m, and the water depth in the southwestern flood land is now usually less than 1 m.After the operation of the PLHP (M2), with the rise of the water level in the lake region, the suitable area for the growth of Vallisneria natans is increased greatly, and the variation is proportional to the DWL during the same period. The increase is most obvious in the low-water-level year (2006), followed by the medium water level year (2018), and becomes insignificant in the high-water-level year (2010). Such a trend is consistent with the previously mentioned variation in the water depth.Since the DWL is relatively small during September and October in 2010, the habitat suitability in this period is changed little between M1 and M2 (Fig. 6). The distribution of the habitat suitability is completely the same between Fig. 6c,f, while the difference in the distributions of habitat suitability between Fig. 6i,l is subtle. As the water level in M1 is relatively low during September and October in 2006 and 2018, the suitable area for the growth of Vallisneria natans in M2 is expanded from northeast to southwest under the influence of the PLHP. In addition, the habitat suitability is increased greatly in Wucheng National Nature Reserve and Nanji National Nature Reserve. Before the operation of the PLHP, there are large areas of dry land in the two reserves, and thus the habitat suitability in these dry areas is 0. After the operation of the PLHP, according the previous research, more than 1/3 of area in the two reserves sees apparent increases in water depth. The water depth is increased to 1 ~ 2 m and thus the habitat suitability is increased to 1.0 in the northeastern part of the reserves. In the southwestern part of the reserves, although the increase of water depth is not significant, most of the lake bed has changed its status from being dry to being wet and the habitat suitability is correspondingly increased from 0 to 0.1 ~ 0.2, as shown in Figs. 6d,e,j,k. This means that the two nature reserves will be more suitable for the growth of Vallisneria natans after the operation of the PLHP, and there it will be easier for Siberian Crane to find food here.Changes in habitat area of Vallisneria natansOn the basis of formula (6), we calculate the habitat area from 1st of September to 31st of October in three typical years, as shown in Fig. 7 and Table 3. In 2006, the habitat area is increased greatly, and the impact of the PLHP on habitat area is most evident in October. Compared with M1, the mean monthly habitat area in M2 is increased by 190.92%. Especially around the 10th of October, the increase can reach about 867.79 km2, accounting for about 1/4 of the total area of the lake. Compared with 2006, the increase of habitat area in 2018 is smaller. The mean monthly habitat area in M2 is increased by 57.07% in September and by 145.27% in October. The largest change occurs in late September, with an increase of about 841.43 km2, which is basically the same as in 2006. While in 2010, compared with M1, the habitat area in M2 is completely the same in September and the average increase in October is only 18.07%, indicating that when the water level is relatively high the operation of the PLHP will make little change to the habitat area.Figure 7Variations of habitat areas of Vallisneria natans in M1 (without the PLHP) and M2 (with the PLHP) and the resulting differences (grey columns).Full size imageTable 3 Mean monthly habitat areas of Vallisneria natans and changes between M1 and M2.Full size tableAfter the operation of the PLHP, the habitat area can reach more than 1000 km2 during the most time of September and October in three typical years, accounting for about 1/3 of the total area of the lake region. In other words, the latest official regulatory scheme of the PLHP is beneficial for the growth of Vallisneria natans in September and October. It means that, whether it is a low-water-level year, medium-water-level year or high-water-level year, before Siberian Crane fly to Poyang Lake for winter, Vallisneria natans will occupy large areas of the flooded land under the regulation of the PLHP. It will ensure that Siberian Crane can consume abundant tubers of Vallisneria natans as food in winter, allowing them to better survive and reproduce in the wetlands of Poyang Lake.In this research, the habitat suitability model of Vallisneria natans in Poyang Lake is established based on the previous research of Chen et al.16. This model only considers the effect of water depth, which mainly influences the growth of aquatic plant by changing the degree of light attenuation34. In fact, temperature and flow velocity can also influence the growth of Vallisneria natans according the relevant studies. However, temperature only plays a decisive role in the germination period38 and is rarely considered as an influencing factor during the growth period of Vallisneria natans. While high flow velocity may adversely affect the growth of Vallisneria natans during the seedling period, adult Vallisneria plants have an extensive root–rhizome system and long ribbon-like leaves to prevent them from being torn apart in rapid water39. Vallisneria natans often begin to sprout in March, reaching the tillering stage in June or July in Poyang Lake region. As a consequence, the temperature and flow velocity will have little effect on the growth of Vallisneria natans during the study period in this research (September and October). Therefore, the water depth can be regarded as the main factor affecting the suitability of Vallisneria natans during the mature period.There have been several scholars who have studied the effect of water depth on the growth of Vallisneria natans in other areas. Xiao et al.40 reported that Vallisneria natans grow rapidly with depths of 110–160 cm, while the growth is severely retarded with a depth of 250 cm. Cao et al.34 carried out experiments in turbid water and reported that the water depth of about 130 cm is most suitable for the growth of Vallisneria natans, while higher water depth will be less favorable for the growth. The inconsistency between these findings and the habitat suitability curve proposed by Chen et al.16 may be explained by the climate differences in different regions, as well as the different degrees of turbidity in water. According to the above analysis, the habitat suitability model of Vallisneria natans in this paper is established based on the habitat suitability curve proposed by Chen et al. (2020) as it represents to the growth characteristics of Vallisneria natans in Poyang Lake region.According to the above analysis, the latest regulatory scheme of the PLHP will effectively increase the water level in the lake region and expand the habitat area of Vallisneria natans, especially in low-water-level years. This finding is different from some of the previous research. Zhu et al.15 used the average water depth during the growing period of Vallisneria natans (from March to October) to reflect the availability of this food resource for Siberian Crane. They found that the PLHP had few influences on Vallisneria natans. This was because the PLHP remain completely open from April to August and it takes effect only in March, September and October during the growing period of Vallisneria natans. Therefore, the average water depth during March to October differs little whether with or without the PLHP, which would certainly underestimate the impact of the PLHP. The present study focuses on September and October, and uses the mean monthly water depth to reflect the habitat suitability of Vallisneria natans. In this period, the natural water level is relatively low in M1, and the operation of the PLHP increases the water level and inundates large areas of otherwise dry land. As a result, the habitat areas of Vallisneria natans observe an increase. More

(ECHA), E. C. A. Know more about the effects of the chemicals we use in Europe (ECHA/PR/16/01). https://echa.europa.eu/de/-/know-more-about-the-effects-of-the-chemicals-we-use-in-europe (2016).Liu, W. J., Nie, H. X., Liang, D. P., Bai, Y. & Liu, H. W. Phospholipid imaging of zebrafish exposed to fipronil using atmospheric pressure matrix-assisted laser desorption ionization mass spectrometry. Talanta https://doi.org/10.1016/j.talanta.2019.120357 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Sparvero, L. J. et al. Mapping of phospholipids by MALDI imaging (MALDI-MSI): Realities and expectations. Chem. Phys. Lipid. 165, 545–562. https://doi.org/10.1016/j.chemphyslip.2012.06.001 (2012).CAS 
Article 

Google Scholar 
Koizumi, S. et al. Imaging mass spectrometry revealed the production of lyso-phosphatidylcholine in the injured ischemic rat brain. Neuroscience 168(1), 219–225. https://doi.org/10.1016/j.neuroscience.2010.03.056 (2010).CAS 
Article 
PubMed 

Google Scholar 
Hankin, J. A. et al. MALDI mass spectrometric imaging of lipids in rat brain injury models. J. Am. Soc. Mass Spectrom. 22(6), 1014–1021. https://doi.org/10.1007/s13361-011-0122-z (2011).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Zhao, C. et al. MALDI-MS imaging reveals asymmetric spatial distribution of lipid metabolites from bisphenol s-induced nephrotoxicity. Anal. Chem. 90(5), 3196–3204. https://doi.org/10.1021/acs.analchem.7b04540 (2018).CAS 
Article 
PubMed 

Google Scholar 
Barbacci, D. C. et al. Mass spectrometric imaging of ceramide biomarkers tracks therapeutic response in traumatic brain injury. ACS Chem. Neurosci. 8(10), 2266–2274. https://doi.org/10.1021/acschemneuro.7b00189 (2017).CAS 
Article 
PubMed 

Google Scholar 
Rompp, A. et al. Histology by mass spectrometry: Label-free tissue characterization obtained from high-accuracy bioanalytical imaging. Angew. Chem. Int. Ed. 49, 3834–3838. https://doi.org/10.1002/anie.200905559 (2010).CAS 
Article 

Google Scholar 
Zemski Berry, K. A. et al. MALDI imaging of lipid biochemistry in tissues by mass spectrometry. Chem. Rev. 111, 6491–6512. https://doi.org/10.1021/cr200280p (2011).CAS 
Article 

Google Scholar 
Cornett, D. S., Reyzer, M. L., Chaurand, P. & Caprioli, R. M. MALDI imaging mass spectrometry: Molecular snapshots of biochemical systems. Nat. Methods 4, 828–833. https://doi.org/10.1038/nmeth1094 (2007).CAS 
Article 
PubMed 

Google Scholar 
Römpp, A. & Spengler, B. Mass spectrometry imaging with high resolution in mass and space. Histochem. Cell Biol. 139, 759–783. https://doi.org/10.1007/s00418-013-1097-6 (2013).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Monroe, E. B. et al. SIMS and MALDI MS imaging of the spinal cord. Proteomics 8(18), 3746-3754. https://doi.org/10.1002/pmic.200800127 (2008).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Chaurand, P., Cornett, D. S., Angel, P. M. & Caprioli, R. M. From whole-body sections down to cellular level, multiscale imaging of phospholipids by MALDI mass spectrometry. Mol. Cell. Proteom. https://doi.org/10.1074/mcp.O110.004259 (2011).Article 

Google Scholar 
Lee, H.-B. & Peart, T. E. Determination of bisphenol A in sewage effluent and sludge by solid-phase and supercritical fluid extraction and gas chromatography/mass spectrometry. J. AOAC Int. 83, 290–298. https://doi.org/10.1093/jaoac/83.2.290 (2000).CAS 
Article 
PubMed 

Google Scholar 
Desbenoit, N., Walch, A., Spengler, B., Brunelle, A. & Römpp, A. Correlative mass spectrometry imaging, applying time-of-flight secondary ion mass spectrometry and atmospheric pressure matrix-assisted laser desorption/ionization to a single tissue section. Rapid Commun. Mass Spectrometry 32, 159–166. https://doi.org/10.1002/rcm.8022 (2018).ADS 
CAS 
Article 

Google Scholar 
Meding, S. et al. Tumor classification of six common cancer types based on proteomic profiling by MALDI imaging. J. Proteome Res. 11, 1996–2003. https://doi.org/10.1021/pr200784p (2012).CAS 
Article 
PubMed 

Google Scholar 
Ritschar, S. et al. Classification of target tissues of Eisenia fetida using sequential multimodal chemical analysis and machine learning. Histochem. Cell Biol. https://doi.org/10.1007/s00418-021-02037-1 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Altshuler, I. et al. An integrated multi-disciplinary approach for studying multiple stressors in freshwater ecosystems: Daphnia as a model organism. Integr. Comp. Biol. 51(4), 623–633. https://doi.org/10.1093/icb/icr103 (2011).CAS 
Article 
PubMed 

Google Scholar 
Bambino, K. & Chu, J. in Zebrafish at the Interface of Development and Disease Research Vol. 124 Current Topics in Developmental Biology (ed K. C. Sadler) 331–367 (2017).Seda, J. & Petrusek, A. Daphnia as a model organism in limnology and aquatic biology: Introductory remarks. J. Limnol. 70, 337–344. https://doi.org/10.4081/jlimnol.2011.337 (2011).Article 

Google Scholar 
de Souza Anselmo, C., Sardela, V. F., de Sousa, V. P. & Pereira, H. M. G. Zebrafish (Danio rerio): A valuable tool for predicting the metabolism of xenobiotics in humans? Comp. Biochem. Physiol. Part C: Toxicol. Pharmacol. 212, 34–46. https://doi.org/10.1016/j.cbpc.2018.06.005 (2018).CAS 
Article 

Google Scholar 
Panula, P. et al. The comparative neuroanatomy and neurochemistry of zebrafish CNS systems of relevance to human neuropsychiatric diseases. Neurobiol. Dis. 40, 46–57. https://doi.org/10.1016/j.nbd.2010.05.010 (2010).CAS 
Article 
PubMed 

Google Scholar 
Korn, H. & Faber, D. S. The Mauthner cell half a century later: A neurobiological model for decision-making?. Neuron 47, 13–28. https://doi.org/10.1016/j.neuron.2005.05.019 (2005).CAS 
Article 
PubMed 

Google Scholar 
Schirmer, E., Schuster, S. & Machnik, P. Bisphenols exert detrimental effects on neuronal signaling in mature vertebrate brains. Commun. Biol. https://doi.org/10.1038/s42003-021-01966-w (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Flößner, D. Book review: Cladocera: The genus Daphnia (including Daphniopsis). Int. Rev. Hydrobiol. 90, 637. https://doi.org/10.1002/iroh.200590003 (2005).Article 

Google Scholar 
OECD. Test No. 211: Daphnia magna Reproduction Test. (2012).Muyssen, B. T. A. & Janssen, C. R. Multigeneration zinc acclimation and tolerance in Daphnia magna: Implications for water-quality guidelines and ecological risk assessment. Environ. Toxicol. Chem. 20, 2053–2060. https://doi.org/10.1002/etc.5620200926 (2001).CAS 
Article 
PubMed 

Google Scholar 
Blewett, T. A. et al. Sublethal and reproductive effects of acute and chronic exposure to flowback and produced water from hydraulic fracturing on the water flea Daphnia magna. Environ. Sci. Technol. 51, 3032–3039. https://doi.org/10.1021/acs.est.6b05179 (2017).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Yang, J. H., Kim, H. J., Lee, S. M., Kim, B. M. & Seo, Y. R. Cadmium-induced biomarkers discovery and comparative network analysis in Daphnia magna. Mol. Cell. Toxicol. 13, 327–336. https://doi.org/10.1007/s13273-017-0036-3 (2017).CAS 
Article 

Google Scholar 
Ferain, A. et al. Body lipid composition modulates acute cadmium toxicity in Daphnia magna adults and juveniles. Chemosphere 205, 328–338. https://doi.org/10.1016/j.chemosphere.2018.04.091 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Ritschar, S., Narayana, V. K. B., Rabus, M. & Laforsch, C. Uncovering the chemistry behind inducible morphological defences in the crustacean Daphniamagna via micro-Raman spectroscopy. Sci. Rep. 10(1), 22408. https://doi.org/10.1038/s41598-020-79755-4 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Machnik, P., Schirmer, E., Glück, L. & Schuster, S. Recordings in an integrating central neuron provide a quick way for identifying appropriate anaesthetic use in fish. Sci. Rep. 8, 17541. https://doi.org/10.1038/s41598-018-36130-8 (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Luzio, A. et al. Copper induced upregulation of apoptosis related genes in zebrafish (Danio rerio) gill. Aquat. Toxicol. 128, 183–189. https://doi.org/10.1016/j.aquatox.2012.12.018 (2013).CAS 
Article 
PubMed 

Google Scholar 
Macirella, R. & Brunelli, E. Morphofunctional alterations in zebrafish (Danio rerio) gills after exposure to mercury chloride. Int. J. Mol. Sci. https://doi.org/10.3390/ijms18040824 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Mansouri, B. & Johari, S. A. Effects of short-term exposure to sublethal concentrations of silver nanoparticles on histopathology and electron microscope ultrastructure of zebrafish (Danio rerio) gills. IJT 10, 15–20. https://doi.org/10.32598/IJT.10.1.60.4 (2016).CAS 
Article 

Google Scholar 
Perez, C. J., Tata, A., de Campos, M. L., Peng, C. & Ifa, D. R. Monitoring toxic ionic liquids in zebrafish (Danio rerio) with desorption electrospray ionization mass spectrometry imaging (DESI-MSI). J. Am. Soc. Mass Spectrom. 28, 1136–1148. https://doi.org/10.1007/s13361-016-1515-9 (2017).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Stutts, W. L. et al. Methods for cryosectioning and mass spectrometry imaging of whole-body zebrafish. J. Am. Soc. Mass Spectrom. 31, 768–772. https://doi.org/10.1021/jasms.9b00097 (2020).CAS 
Article 
PubMed 

Google Scholar 
Purves, D. & Williams, S. M. Neuroscience. 2nd edition. Vol. Chapter 11, Vision: The Eye (Sinauer Associates, 2001).
Google Scholar 
Strungaru, S. A. et al. Toxicity and chronic effects of deltamethrin exposure on zebrafish (Danio rerio) as a reference model for freshwater fish community. Ecotoxicol. Environ. Saf. 171, 854–862. https://doi.org/10.1016/j.ecoenv.2019.01.057 (2019).CAS 
Article 
PubMed 

Google Scholar 
Mishra, A. & Devi, Y. Histopathological alterations in the brain (optic tectum) of the fresh water teleost Channa punctatus in response to acute and subchronic exposure to the pesticide Chlorpyrifos. Acta Histochem. 116, 176–181. https://doi.org/10.1016/j.acthis.2013.07.001 (2014).CAS 
Article 
PubMed 

Google Scholar 
Jia, W., Mao, L., Zhang, L., Zhang, Y. & Jiang, H. Effects of two strobilurins (azoxystrobin and picoxystrobin) on embryonic development and enzyme activities in juveniles and adult fish livers of zebrafish (Danio rerio). Chemosphere 207, 573–580. https://doi.org/10.1016/j.chemosphere.2018.05.138 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Seyoum, A., Pradhan, A., Jass, J. & Olsson, P. E. Perfluorinated alkyl substances impede growth, reproduction, lipid metabolism and lifespan in Daphnia magna. Sci. Total Environ. https://doi.org/10.1016/j.scitotenv.2020.139682 (2020).Article 
PubMed 

Google Scholar 
Scanlan, L. D. et al. Gene transcription, metabolite and lipid profiling in eco-indicator Daphnia magna indicate diverse mechanisms of toxicity by legacy and emerging flame-retardants. Environ. Sci. Technol. 49, 7400–7410. https://doi.org/10.1021/acs.est.5b00977 (2015).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Heinlaan, M. et al. Changes in the Daphnia magna midgut upon ingestion of copper oxide nanoparticles: A transmission electron microscopy study. Water Res. 45, 179–190. https://doi.org/10.1016/j.watres.2010.08.026 (2011).CAS 
Article 
PubMed 

Google Scholar 
Abe, T., Saito, H., Niikura, Y., Shigeoka, T. & Nakano, Y. Embryonic development assay with Daphnia magna: Application to toxicity of aniline derivatives. Chemosphere 45, 487–495. https://doi.org/10.1016/s0045-6535(01)00049-2 (2001).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Sengupta, N., Gerard, P. D. & Baldwin, W. S. Perturbations in polar lipids, starvation survival and reproduction following exposure to unsaturated fatty acids or environmental toxicants in Daphnia magna. Chemosphere 144, 2302–2311. https://doi.org/10.1016/j.chemosphere.2015.11.015 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Huber, K. et al. Approaching cellular resolution and reliable identification in mass spectrometry imaging of tryptic peptides. Anal. Bioanal. Chem. 410, 5825–5837. https://doi.org/10.1007/s00216-018-1199-z (2018).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
White, R. M. et al. Transparent adult zebrafish as a tool for in vivo transplantation analysis. Cell Stem Cell 2, 183–189. https://doi.org/10.1016/j.stem.2007.11.002 (2008).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Nagayoshi, S. et al. Insertional mutagenesis by the Tol2 transposon-mediated enhancer trap approach generated mutations in two developmental genes: tcf7 and synembryn-like. Development 135, 159–169. https://doi.org/10.1242/dev.009050 (2008).CAS 
Article 
PubMed 

Google Scholar 
Perciedu Sert, N. et al. The ARRIVE guidelines 2.0: Updated guidelines for reporting animal research. Exp. Physiol. 105, 1459–1466. https://doi.org/10.1113/EP088870 (2020).Article 

Google Scholar 
Elendt, B. P. Selenium deficiency in Crustacea. Protoplasma 154, 25–33. https://doi.org/10.1007/BF01349532 (1990).CAS 
Article 

Google Scholar 
Sud, M. et al. LMSD: LIPID MAPS structure database. Nucleic Acids Res. 35, D527–D532. https://doi.org/10.1093/nar/gkl838 (2007).CAS 
Article 
PubMed 

Google Scholar 
Race, A. M., Styles, I. B. & Bunch, J. Inclusive sharing of mass spectrometry imaging data requires a converter for all. J. Proteom. 75, 5111–5112. https://doi.org/10.1016/j.jprot.2012.05.035 (2012).CAS 
Article 

Google Scholar 
Robichaud, G., Garrard, K. P., Barry, J. A. & Muddiman, D. C. MSiReader: An open-source interface to view and analyze high resolving power MS imaging files on Matlab platform. J. Am. Soc. Mass Spectrom. 24, 718–721. https://doi.org/10.1007/s13361-013-0607-z (2013).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar  More

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

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

Dordevic, B. & Dasic, T. Water storage reservoirs and their role in the development, utilization and protection of catchment. SPATIUM Int. Rev. 24, 9–15 (2011).
Google Scholar 
European Community. Directive 2000/60/EC of the European Parliament and of the Council of 23 October 2000 establishing a framework for Community action in the field of water policy. Official J. Eur. Union. 5, L327 (2000).

Google Scholar 
Riley, W. D. et al. Small water bodies in Great Britain and Ireland: Ecosystem function, human-generated degradation, and options for restorative action. Sci. Total Environ. 645, 1598–1616. https://doi.org/10.1016/j.scitotenv.2018.07.243 (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Kujawa, K., Arczyńska-Chudy, E., Janku, K. & Mana, M. Effect of buffer zone structure on diversity of aquatic vegetation in farmland water bodies. Pol. J. Ecol. 68(4), 263–282. https://doi.org/10.3161/15052249PJE2020.68.4.001 (2021).Article 

Google Scholar 
Akasaka, M., Takamura, N., Mitsuhashi, H. & Kadono, Y. Effects of land use on aquatic macrophyte diversity and water quality of ponds. Fresh Biol. 55, 909–922 (2010).CAS 
Article 

Google Scholar 
Pukacz, A., Pełechaty, M., Pełechata, A. & Siepak, M. The differential cover of submerged vegetation vs habitat conditions in the lakes of the Lubuskie Region. Limnol. Rev. 7(2), 95–100 (2007).
Google Scholar 
Scheffer, M., Hosper, S. H., Meijer, M. L., Moss, B. & Jeppesen, E. Alternative equilibria in shallow lakes. Trends Ecol. Evol. 8, 275–279 (1993).CAS 
Article 

Google Scholar 
Scheffer, M. & Jeppesen, E. Alternative stable states. Ecol. Stud. 131, 397–405 (1998).Article 

Google Scholar 
Beck, M. W., Tomcko, C., Valley, R. D. & Staples, D. F. Analysis of macrophyte indicator variation as a function of sampling, temporal and stressor effects. Ecol. Indic. 46, 322–335 (2014).Article 

Google Scholar 
Celewicz-Gołdyn, S. & Kuczyńska-Kippen, N. Ecological value of macrophyte cover in creating habitat for microalgae (diatoms) and zooplankton (rotifers and crustaceans) in small field and forest water bodies. PLoS ONE 12(5), e0177317. https://doi.org/10.1371/journal.pone.0177317 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Wilk-Woźniak, E. et al. Effects of the environs of waterbodies on aquatic plants in oxbow lakes (habitat 3150). Ecol. Ind. 98, 736–742. https://doi.org/10.1016/j.ecolind.2018.11.025 (2019).CAS 
Article 

Google Scholar 
Bedla, D. & Misztal, A. Changeability of chemistry of small water reservoirs with diversified use structure of the adjoining areas. Annu. Set Environ. Protect. 16, 421–439 (2014) (ISSSN 1506-218X (in Polish)).
Google Scholar 
Mozgawa, J. Photointerpretation analysis of landscape structure in lake watersheds of Suwałki Landscape. Park. Ekol. Pol. 41, 53–74 (1993).
Google Scholar 
Mioduszewski, W. Small water reservoirs on agricultural areas. Wieś Jutra Zakład Zasobów Wodnych Instytut Melioracji i Użytków Zielonych Falenty. 10(123), 32–34 (2008) (in Polish).
Google Scholar 
Gołdyn, B., Kowalczewska-Madura, K. & Celewicz-Gołdyn, S. Drought and deluge: Influence of environmental factors on water quality of kettle holes in two subsequent years with different precipitation. Limnologica 54, 14–22. https://doi.org/10.1016/j.limno.2015.07.002 (2015).CAS 
Article 

Google Scholar 
Bylak, A. et al. Small stream catchment in a developing city context: The importance of land cover changes on the ecological status of streams and the possibilities for providing ecosystem services. Sci. Total Environ. 815, 151974. https://doi.org/10.1016/j.scitotenv.2021.151974 (2022).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Goyal, V. C. et al. Ecological health and water quality of village ponds in the subtropics limiting their use for water supply and groundwater recharge. J. Environ. Manage. 277, 111450. https://doi.org/10.1016/j.jenvman.2020.111450 (2021).CAS 
Article 
PubMed 

Google Scholar 
Kuczyńska-Kippen, N., Spoljar, M., Mleczek, M. & Zhang, Ch. Elodeides, but not helophytes, increase community diversity and reduce trophic state: Case study with rotifer indices in field ponds. Ecol. Ind. 128, 107829. https://doi.org/10.1016/j.ecolind.2021.107829 (2021).Article 

Google Scholar 
Novikmec, M. et al. Ponds and their catchments: size relationships and influence of land use across multiple spatial scales. Hydrobiologia 774, 155–166. https://doi.org/10.1007/s10750-015-2514-8 (2016).Article 

Google Scholar 
Dudzińska, A., Szpakowska, B. & Pajchrowska, M. Influence of land development on the ecological status of small water bodies. Oceanol. Hydrobiol. Stud. 49(4), 345–353. https://doi.org/10.1515/ohs-2020-0030 (2020).Article 

Google Scholar 
Szpakowska, B. Occurrence and Role of Organic Compounds Dissolved in Surface and Ground Waters of Agricultural Landscape (Wydawnictwo Uniwersytetu Mikołaja Kopernika, 1999) (in Polish).
Google Scholar 
Wysocka-Czubaszek, A. & Banaszuk, P. Migration of nitrogen compounds and the riparian zones in the Upper Narew Valley. Acta Agroph. 2(1), 349–354 (2003) (in Polish).
Google Scholar 
Szpakowska, B., Świerk, D., Pajchrowska, M. & Gołdyn, R. Verifying the usefulness of macrophytes as an indicator of the status of small waterbodies. Sci. Total Environ. 798, 149279. https://doi.org/10.1016/j.scitotenv.2021.149279 (2021).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Hermanowicz, W., Dojlido, J., Dożańska, W., Koziorowski, B. & Zerbe, J. Physico-Chemical Examinations of Water and Sludge (Arkady, 1999) (in Polish).
Google Scholar 
Aynur, M., Liming, N. & Fang, Y. Spatial evaluation of environmental suitability for human settlement of Kashgar, Northwest China. Fresenius Environ. Bull. 27(9), 5899–5907 (2018).
Google Scholar 
Pham, L., Brabyn, L. & Ashrof, S. Combining Quick Bird, LiDAR and GIS topography indices to identify a single native tree species in a complex landscape using an object-based classification approach. Int. J. Appl. Earth Obs. Geoinf. 50, 187–197. https://doi.org/10.1016/j.jag.2016.03.015 (2016).ADS 
Article 

Google Scholar 
Lepš, J. & Šmilauer, P. Multivariate Analysis of Ecological Data Using CANOCO (Cambridge Univeristy Press, 2003) (ISBN 9780521891080).Book 

Google Scholar 
Mohamed, Z. A. Macrophytes-cyanobacteria allelopathic interactions and their implications for water resources management—A review. Limnologica 63, 122–132. https://doi.org/10.1016/j.limno.2017.02.006 (2017).CAS 
Article 

Google Scholar 
Koc, J. Problems of small reservoirs protection on rural areas. In Problemy ochrony i użytkowania obszarów wiejskich o dużych walorach przyrodniczych (eds Radwan, S. & Lorkiewicz, Z.) 151–156 (UMCS, 2000).
Google Scholar 
Bell, V. A. et al. Long term simulations of macronutrients (C, N and P) in UK freshwaters. Sci. Total Environ. 776, 145813. https://doi.org/10.1016/j.scitotenv.2021.145813 (2021).ADS 
CAS 
Article 

Google Scholar 
Winton, R. S. et al. Anthropogenic influences on Zambian water quality: Hydropower and land-use change. Environ. Sci. Process. Impacts. 23, 981–994. https://doi.org/10.1039/D1EM00006C (2021).CAS 
Article 
PubMed 

Google Scholar 
Wang, Y. et al. Effect and risk assessment of animal manure pollution on Huaihe River Basin, China. Chin. Geogr. Sci. 31, 751–764. https://doi.org/10.1007/s11769-021-1222-8 (2021).Article 

Google Scholar 
Lawniczak-Malińska, A., Ptak, M., Celewicz, S. & Choiński, A. Impact of lake morphology and shallowing on the rate of overgrowth in hard-water eutrophic lakes. Water 10(12), 1827. https://doi.org/10.3390/w10121827 (2018).CAS 
Article 

Google Scholar 
Bosiacka, B., Pacewicz, K. & Pieńkowski, P. Spatial analysis of plant species distribution among small water bodies in an agricultural landscape. Acta Agrobot. 61(2), 93–101. https://doi.org/10.5586/aa.2008.037 (2008).Article 

Google Scholar 
Joniak, T., Kuczyńska-Kippen, N. & Gąbka, M. Effect of agricultural landscape characteristics on the hydrobiota structure in small water bodies. Hydrobiologia 793, 121–133. https://doi.org/10.1007/s10750-016-2913-5 (2017).CAS 
Article 

Google Scholar 
Barling, R. D. & Moore, I. D. Role of buffer strips in management of waterway pollution: A review. Environ. Manage. 18, 543–558. https://doi.org/10.1007/BF02400858 (1994).ADS 
Article 

Google Scholar 
Blanco-Cangui, H., & Lal, R. Buffer strips. In Principles of Soil Conservation and Management (eds Blanco, H. & Lal, R.) 223–257 (Springer, Netherlands, 2008).
Google Scholar 
Borin, M., Passoni, M., Thiene, M. & Tempesta, T. Multiple functions of buffer strips in farming areas. Eur. J. Agron. 32, 103–111. https://doi.org/10.1016/j.eja.2009.05.003 (2010).Article 

Google Scholar 
Ghobrial, M. G., Nassr, H. S. & Kamil, A. W. Bioactivity effect of two macrophyte extracts on growth performance of two bloom-forming cyanophytes. Egypt. J. Aquat. Res. 41(1), 69–81. https://doi.org/10.1016/j.ejar.2015.01.001 (2015).Article 

Google Scholar 
Lacas, J. G., Voltz, M., Gouy, V., Carluer, N. & Gril, J. J. Using grassed strips to limit pesticide transfer to surface water: A review. Agron. Sustain. Dev. 25, 253–266 (2005).CAS 
Article 

Google Scholar 
Świerk, D. & Szpakowska, B. An ecosystem valuation method for small water bodies. Ecol. Chem. Eng. S. 20(2), 397–418. https://doi.org/10.2478/eces-2013-0029 (2013).Article 

Google Scholar 
Marszałek, M., Kowalski, Z. & Makara, A. The possibility of contamination of water-soil environment as a results of the use of pig slurry. Ecol. Chem. Eng. S. 26(2), 313–330. https://doi.org/10.1515/eces-219-0022 (2019).Article 

Google Scholar 
Micek, G., Górecki, J. & Neo, H. Relations: Company and its local milieu in the context of foreign direct investment in Polish pig industry. In Człowiek i Środowisko (eds Górka, Z. & Zborowski, A.) 297–308 (UJ Kraków, 2009).
Google Scholar 
OECD. Agriculture Trade and The Environment: The Pig Sector (OECD, 2003).
Google Scholar 
Barałkiewicz, D. et al. Storm water contamination and its effect on the quality of urban surface waters. Environ. Monit. Assess. 186(10), 6789–6803. https://doi.org/10.1007/s10661-014-3889-0 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
So Fijanic, A., Hulley, M. & Loock, D. Stormwater quality assessment and management for the town of jasper in Alberta, Canada. Water Qual. Res. J. Can. 56(1), 45–56. https://doi.org/10.2166/wqrj.2021.012 (2021).CAS 
Article 

Google Scholar 
Zubala, T. Assessment of the variability of rainwater quality and the functioning of retention reservoirs in the urban area. Rocznik Ochrona Środowiska 22(2), 840–858 (2020).
Google Scholar 
Czerniawski, R., Sługocki, L., Krepski, T., Wilczak, A. & Pietrzak, K. Spatial changes in invertebrate structures as a factor of strong human activity in the bed and catchment area of a small urban stream. Water 12(3), 913. https://doi.org/10.3390/w12030913 (2020).CAS 
Article 

Google Scholar 
Gołdyn, R. et al. Influence of stormwater runoff on macroinvertebrates in a small urban river and a reservoir. Sci. Total Environ. 625, 743–751. https://doi.org/10.1016/j.scitotenv.2017.12.324 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Cao, Y., Zhi, Y., Jeppesen, E. & Li, W. Species-specific responses of submerged macrophytes to simulated extreme precipitation: A mesocosm study. Water 11(6), 1160. https://doi.org/10.3390/w11061160 (2019).CAS 
Article 

Google Scholar 
Hilt, S. et al. Restoration of submerged vegetation in shallow eutrophic lakes—A guideline and state of the art in Germany. Limnologica 36, 155–171 (2006).CAS 
Article 

Google Scholar 
Ryszkowski, L. & Kędziora, A. Modification of water flows and nitrogen fluxes by shelterbelts. Ecol. Eng. 29(4), 388–400. https://doi.org/10.1016/j.ecoleng.2006.09.023 (2007).Article 

Google Scholar  More

Schimel, J. P., Bennett, J. & Fierer, N. Microbial community composition and soil nitrogen cycling: is there really a connection? In Biological Diversity and Function in Soils Ecological Reviews (eds Bardgett, R. et al.) 171–188 (Cambridge University Press, 2005).
Google Scholar 
Hatzenpichler, R. Diversity, physiology, and niche differentiation of ammonia-oxidizing archaea. Appl. Environ. Microbiol. 78, 7501–7510 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
Kowalchuk, G. A. & Stephen, J. R. Ammonia-oxidizing bacteria: A model for molecular microbial ecology. Annu. Rev. Microbiol. 55, 485–529 (2001).CAS 
PubMed 

Google Scholar 
Daims, H. et al. Complete nitrification by Nitrospira bacteria. Nature 528, 504–509 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
Shi, X. et al. Niche separation of comammox Nitrospira and canonical ammonia oxidizers in an acidic subtropical forest soil under long-term nitrogen deposition. Soil Biol. Biochem. 126, 114–122 (2018).CAS 

Google Scholar 
van Kessel, M. A. H. J. et al. Complete nitrification by a single microorganism. Nature 528, 555–559 (2015).PubMed 
PubMed Central 

Google Scholar 
Wang, X. et al. Comammox bacterial abundance, activity, and contribution in agricultural rhizosphere soils. Sci. Total Environ. 727, 138563 (2020).CAS 
PubMed 

Google Scholar 
Wang, F., Liang, X., Ma, S., Liu, L. & Wang, J. Ammonia-oxidizing archaea are dominant over comammox in soil nitrification under long-term nitrogen fertilization. J. Soils Sediments 21, 1800–1814 (2021).CAS 

Google Scholar 
Rotthauwe, J.-H., Witzel, K.-P. & Liesack, W. The ammonia monooxygenase structural gene amoA as a functional marker: Molecular fine-scale analysis of natural ammonia-oxidizing populations. Appl. Environ. Microbiol. 63, 4704–4712 (1997).CAS 
PubMed 
PubMed Central 

Google Scholar 
Spang, A. et al. Distinct gene set in two different lineages of ammonia-oxidizing archaea supports the phylum Thaumarchaeota. Trends Microbiol. 18, 331–340 (2010).CAS 
PubMed 

Google Scholar 
Schleper, C. & Nicol, G. W. Ammonia-oxidising archaea—Physiology, ecology and evolution. Adv. Microb. Physiol. 57, 1–41 (2010).CAS 
PubMed 

Google Scholar 
Prosser, J. I. & Nicol, G. W. Archaeal and bacterial ammonia-oxidisers in soil: The quest for niche specialisation and differentiation. Trends Microbiol. 20, 523–531 (2012).CAS 
PubMed 

Google Scholar 
Amin, S. A. et al. Copper requirements of the ammonia-oxidizing archaeon Nitrosopumilus maritimus SCM1 and implications for nitrification in the marine environment. Limnol. Oceanogr. 58, 2037–2045 (2013).CAS 

Google Scholar 
Jenkins, S. N., Murphy, D. V., Waite, I. S., Rushton, S. P. & O’Donnell, A. G. Ancient landscapes and the relationship with microbial nitrification. Sci. Rep. 6, 30733 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Gubry-Rangin, C. et al. Niche specialization of terrestrial archaeal ammonia oxidizers. Proc. Natl. Acad. Sci. U.S.A. 108, 21206–21211 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
Lehtovirta-Morley, L. E., Stoecker, K., Vilcinskas, A., Prosser, J. I. & Nicol, G. W. Cultivation of an obligate acidophilic ammonia oxidizer from a nitrifying acid soil. Proc. Natl. Acad. Sci. U.S.A. 108, 15892–15897 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
Banning, N. C., Maccarone, L. D., Fisk, L. M. & Murphy, D. V. Ammonia-oxidising bacteria not archaea dominate nitrification activity in semi-arid agricultural soil. Sci. Rep. 5, 11146 (2015).PubMed 
PubMed Central 

Google Scholar 
Di, H. J. et al. Ammonia-oxidizing bacteria and archaea grow under contrasting soil nitrogen conditions. FEMS Microbiol. Ecol. 72, 386–394 (2010).CAS 
PubMed 

Google Scholar 
Wang, J., Wang, J., Rhodes, G., He, J. Z. & Ge, Y. Adaptive responses of comammox Nitrospira and canonical ammonia oxidizers to long-term fertilizations: Implications for the relative contributions of different ammonia oxidizers to soil nitrogen cycling. Sci. Total Environ. 668, 224–233 (2019).CAS 
PubMed 

Google Scholar 
Ouyang, Y., Evans, S. E., Friesen, M. L. & Tiemann, L. K. Effect of nitrogen fertilization on the abundance of nitrogen cycling genes in agricultural soils: A meta-analysis of field studies. Soil Biol. Biochem. 127, 71–78 (2018).CAS 

Google Scholar 
Verhamme, D. T., Prosser, J. I. & Nicol, G. W. Ammonia concentration determines differential growth of ammonia-oxidising archaea and bacteria in soil microcosms. ISME J. 5, 1067–1071 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
Wardle, D. A. Controls of temporal variability of the soil microbial biomass: A global-scale synthesis. Soil Biol. Biochem. 30, 1627–1637 (1998).CAS 

Google Scholar 
Adair, K. L. & Schwartz, E. Evidence that ammonia-oxidizing archaea are more abundant than ammonia-oxidizing bacteria in semiarid soils of northern Arizona, USA. Microb. Ecol. 56, 420–426 (2008).CAS 
PubMed 

Google Scholar 
Taylor, A. E., Zeglin, L. H., Wanzek, T. A., Myrold, D. D. & Bottomley, P. J. Dynamics of ammonia-oxidizing archaea and bacteria populations and contributions to soil nitrification potentials. ISME J. 6, 2024–2032 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
Hayatsu, M., Katsuyama, C. & Tago, K. Overview of recent researches on nitrifying microorganisms in soil. Soil Sci. Plant Nutr. 67, 1–14 (2021).
Google Scholar 
Sher, Y., Zaady, E. & Nejidat, A. Spatial and temporal diversity and abundance of ammonia oxidizers in semi-arid and arid soils: Indications for a differential seasonal effect on archaeal and bacterial ammonia oxidizers. FEMS Microbiol. Ecol 86, 544–556 (2013).CAS 
PubMed 

Google Scholar 
Stopnišek, N. et al. Thaumarchaeal ammonia oxidation in an acidic forest peat soil is not influenced by ammonium amendment. Appl. Environ. Microbiol. 76, 7626–7634 (2010).PubMed 
PubMed Central 

Google Scholar 
Habteselassie, M. Y., Xu, L. & Norton, J. M. Ammonia-oxidizer communities in an agricultural soil treated with contrasting nitrogen sources. Front. Microbiol. 4, 326 (2013).PubMed 
PubMed Central 

Google Scholar 
Wang, C. et al. Climate change amplifies gross nitrogen turnover in montane grasslands of Central Europe in both summer and winter seasons. Glob. Change Biol. 22, 2963–2978 (2016).
Google Scholar 
Wessén, E., Nyberg, K., Jansson, J. K. & Hallin, S. Responses of bacterial and archaeal ammonia oxidizers to soil organic and fertilizer amendments under long-term management. Appl. Soil Ecol. 45, 193–200 (2010).
Google Scholar 
Kong, A. Y. Y., Hristova, K., Scow, K. M. & Six, J. Impacts of different N management regimes on nitrifier and denitrifier communities and N cycling in soil microenvironments. Soil Biol. Biochem. 42, 1523–1533 (2010).CAS 
PubMed 
PubMed Central 

Google Scholar 
Harrison, P. & Pearce, F. AAAS Atlas of Population & Environment 204 (University of California Press, 2000).
Google Scholar 
Reynolds, J. F., Maestre, F. T., Kemp, P. R., Smith, D. M. S. & Lambin, E. F. Natural and human dimensions of land degradation in drylands: Causes and consequences. In Terrestrial Ecosystems in a Changing World Global Change—The IGBP Series (eds Canadell, J. G. et al.) 247–258 (Springer, 2007).
Google Scholar 
McArthur, W. M. Reference Soils of South-Western Australia 2nd edn. (Department of Agriculture, 2004).
Google Scholar 
Barton, L., Murphy, D. V. & Butterbach-Bahl, K. Influence of crop rotation and liming on greenhouse gas emissions from a semi-arid soil. Agric. Ecosyst. Environ. 167, 23–32 (2013).CAS 

Google Scholar 
Barton, L., Hoyle, F. C., Stefanova, K. T. & Murphy, D. V. Incorporating organic matter alters soil greenhouse gas emissions and increases grain yield in a semi-arid climate. Agric. Ecosyst. Environ. 231, 320–330 (2016).CAS 

Google Scholar 
Gubry-Rangin, C., Nicol, G. W. & Prosser, J. I. Archaea rather than bacteria control nitrification in two agricultural acidic soils. FEMS Microbiol. Ecol. 74, 566–574 (2010).CAS 
PubMed 

Google Scholar 
Gleeson, D. B. et al. Response of ammonia oxidizing archaea and bacteria to changing water filled pore space. Soil Biol. Biochem. 42, 1888–1891 (2010).CAS 

Google Scholar 
O’Sullivan, C. A., Wakelin, S. A., Fillery, I. R. P. & Roper, M. M. Factors affecting ammonia-oxidising microorganisms and potential nitrification rates in southern Australian agricultural soils. Soil Res. 51, 240–252 (2013).
Google Scholar 
Zhang, L.-M., Hu, H.-W., Shen, J.-P. & He, J.-Z. Ammonia-oxidizing archaea have more important role than ammonia-oxidizing bacteria in ammonia oxidation of strongly acidic soils. ISME J. 6, 1032–1045 (2012).CAS 
PubMed 

Google Scholar 
Wang, F. et al. Responses of soil ammonia-oxidizing bacteria and archaea to short-term warming and nitrogen input in a semi-arid grassland on the Loess Plateau. Eur. J. Soil Biol. 102, 103267 (2021).CAS 

Google Scholar 
Bolland, M. D. A. & Brennan, R. F. Phosphorus, copper and zinc requirements of no-till wheat crops and methods of collecting soil samples for soil testing. Aust. J. Exp. Agric. 46, 1051–1059 (2006).CAS 

Google Scholar 
Gilkes, B., Lee, S. & Singh, B. The imprinting of aridity upon a lateritic landscape: An illustration from southwestern Australia. C. R. Geosci. 335, 1207–1218 (2003).
Google Scholar 
Hoyle, F. C. & Murphy, D. V. Influence of organic residues and soil incorporation on temporal measures of microbial biomass and plant available nitrogen. Plant Soil 347, 53–64 (2011).CAS 

Google Scholar 
Noy-Meir, I. Desert ecosystems: Environment and producers. Annu. Rev. Ecol. Syst. 4, 25–51 (1973).
Google Scholar 
Petersen, D. G. et al. Abundance of microbial genes associated with nitrogen cycling as indices of biogeochemical process rates across a vegetation gradient in Alaska. Environ. Microbiol. 14, 993–1008 (2012).CAS 
PubMed 

Google Scholar 
Fisk, L. M., Barton, L., Jones, D. L., Glanville, H. C. & Murphy, D. V. Root exudate carbon mitigates nitrogen loss in a semi-arid soil. Soil Biol. Biochem. 88, 380–389 (2015).CAS 

Google Scholar 
Murphy, D. V., Sparling, G. P., Fillery, I. R. P., McNeill, A. M. & Braunberger, P. Mineralisation of soil organic nitrogen and microbial respiration after simulated summer rainfall events in an agricultural soil. Aust. J. Soil Res. 36, 231–246 (1998).
Google Scholar 
Anderson, G. C., Fillery, I. R. P., Dunin, F. X., Dolling, P. J. & Asseng, S. Nitrogen and water flows under pasture–wheat and lupin–wheat rotations in deep sands in Western Australia 2. Drainage and nitrate leaching. Aust. J. Agric. Res. 49, 345–361 (1998).CAS 

Google Scholar 
Nicholls, N. Local and remote causes of the southern Australian autumn-winter rainfall decline, 1958–2007. Clim. Dyn. 34, 835–845 (2010).
Google Scholar 
Delworth, T. L. & Zeng, F. Regional rainfall decline in Australia attributed to anthropogenic greenhouse gases and ozone levels. Nat. Geosci. 7, 583 (2014).CAS 

Google Scholar 
Alexander, L. V. et al. Trends in Australia’s climate means and extremes: A global context. Aust. Meteorol. Mag. 56, 1–18 (2007).
Google Scholar 
Austin, A. T. et al. Water pulses and biogeochemical cycles in arid and semiarid ecosystems. Oecologia 141, 221–235 (2004).PubMed 

Google Scholar 
Isbell, R. F. The Australian Soil Classification 2nd edn. (CSIRO Publishing, 2002).
Google Scholar 
IUSS Working Group WRB. World Reference Base for Soil Resources 2006, First Update 2007 203 (FAO, 2007).
Google Scholar 
Brookes, P. C., Landman, A., Pruden, G. & Jenkinson, D. S. Chloroform fumigation and the release of soil nitrogen: A rapid direct extraction method to measure microbial biomass nitrogen in soil. Soil Biol. Biochem. 17, 837–842 (1985).CAS 

Google Scholar 
Wu, J., Joergensen, R. G., Pommerening, B., Chaussod, R. & Brookes, P. C. Measurement of soil microbial biomass C by fumigation-extraction—An automated procedure. Soil Biol. Biochem. 22, 1167–1169 (1990).CAS 

Google Scholar 
Krom, M. D. Spectrophotometric determination of ammonia: A study of a modified Berthelot reaction using salicylate and dichloroisocyanurate. Analyst 105, 305–316 (1980).CAS 

Google Scholar 
Kamphake, L. J., Hannah, S. A. & Cohen, J. M. Automated analysis for nitrate by hydrazine reduction. Water Res. 1, 205–216 (1967).CAS 

Google Scholar 
Keeney, D. R. & Bremner, J. M. Comparison and evaluation of laboratory methods of obtaining an index of soil nitrogen availability. Agron. J. 58, 498–503 (1966).CAS 

Google Scholar 
Waring, S. A. & Bremner, J. M. Ammonium production in soil under waterlogged conditions as an index of nitrogen availability. Nature 201, 951–952 (1964).CAS 

Google Scholar 
Griffiths, R. I., Whiteley, A. S., O’Donnell, A. G. & Bailey, M. J. Rapid method for coextraction of DNA and RNA from natural environments for analysis of ribosomal DNA- and rRNA-based microbial community composition. Appl. Environ. Microbiol. 66, 5488–5491 (2000).CAS 
PubMed 
PubMed Central 

Google Scholar 
Francis, C. A., Roberts, K. J., Beman, J. M., Santoro, A. E. & Oakley, B. B. Ubiquity and diversity of ammonia-oxidizing archaea in water columns and sediments of the ocean. Proc. Natl. Acad. Sci. U.S.A. 102, 14683–14688 (2005).CAS 
PubMed 
PubMed Central 

Google Scholar 
Barton, L., Gleeson, D. B., Maccarone, L. D., Zúñiga, L. P. & Murphy, D. V. Is liming soil a strategy for mitigating nitrous oxide emissions from semi-arid soils? Soil Biol. Biochem. 62, 28–35 (2013).CAS 

Google Scholar 
Akaike, H. Likelihood of a model and information criteria. J. Econom. 16, 3–14 (1981).MATH 

Google Scholar 
Cresswell, H. P. & Hamilton, G. J. Bulk density and pore space relations. In Soil Physical Measurement and Interpretation for Land Evaluation (eds McKenzie, N. et al.) 35–58 (CSIRO Publishing, 2002).
Google Scholar 
Rayment, G. E. & Lyons, D. J. Soil Chemical Methods—Australasia 495 (CSIRO Publishing, 2011).
Google Scholar  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

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