Ecology

1.
Greenfield, H. J. The secondary products revolution: The past, the present and the future. World Archaeol. 42, 29–54 (2010).
Google Scholar 
2.
Marciniak, A. The secondary products revolution: Empirical evidence and its current zooarchaeological critique. J. World Prehist. 24, 117–130 (2011).
Google Scholar 

3.
Sherratt, A. The secondary exploitation of animals in the Old World. World Archaeol. 15, 90–104 (1983).
Google Scholar 

4.
Patel, A. K. & Meadow, R. H. South Asian contribution to animal domestication and pastoralism: Bones, genes and archaeology. In The Oxford Handbook of Zooarchaeology (eds Albarella, U. et al.) 1–27 (Oxford University Press, Oxford, 2017). https://doi.org/10.1093/oxfordhb/9780199686476.001.0001.
Google Scholar 

5.
Meadow, R. H. & Patel, A. K. Prehistoric pastoralism in northwestern South Asia from the Neolithic through the Harappan period. In Indus Ethnobiology: New perspective from the field (eds Weber, S. A. & Belcher, W. R.) 65–94 (Lexingtion Books, Lanham, 2003).
Google Scholar 

6.
Wright, R. P. The Ancient Indus: Urbanism, Economy, and Society (Cambridge University Press, Cambridge, 2010).
Google Scholar 

7.
Fuller, D. Q. Indus and non-Indus agricultural traditions: Local developments and crop adoptions on the Indian peninsula. In Indus Ethnobiology (eds Weber, S. & Belcher, W. R.) 243–396 (Lexington Books, Lanham, 2003).
Google Scholar 

8.
Meadow, R. H. Pre- and Proto-Historic agriculture and pastoral transformations in northwestern and South Asia. Rev. Archaeol. 19, 12–21 (1998).
Google Scholar 

9.
Pokharia, A. K. et al. Archaeobotany and archaeology at Kanmer, a Harappan site in Kachchh, Gujarat: Evidence for adaptation in response to climatic variability. Curr. Sci. 100, 1833–1846 (2011).
Google Scholar 

10.
Rissman, P. C. Migratory Pastoralism in Western India in the Second Millennium B.C.: The Evidence from Oriyo Timbo (Chiroda). (University of Microfilms International, 1985).

11.
Weber, S., Kashyap, A. & Harriman, D. Does size matter: The role and significance of cereal grains in the Indus Civilization. Archaeol. Anthropol. Sci. 2, 35–43 (2010).
Google Scholar 

12.
Weber, S. A. Plants and Harappan subsistence: An example of stability and change from Rojdi (Oxford & IBH and American Institute of Indian Studies, Oxford, 1991).
Google Scholar 

13.
Goyal, P. et al. Subsistence system, paleoecology, and 14C chronology at Kanmer, a Harappan site in Gujarat India. Radiocarbon 55, 141–150 (2013).
CAS  Google Scholar 

14.
Pokharia, A. K., Kharakwal, J. S. & Srivastava, A. Archaeobotanical evidence of millets in the Indian subcontinent with some observations on their role in the Indus Civilization. J. Archaeol. Sci. 42, 442–455 (2014).
Google Scholar 

15.
Pokharia, A. K. et al. Altered cropping pattern and cultural continuation with declined prosperity following abrupt and extreme arid event at ~4,200 yrs BP: Evidence from an Indus archaeological site Khirsara, Gujarat, western India. PLoS ONE 12, 1–17 (2017).
Google Scholar 

16.
Chase, B. Social change at the Harappan settlement of Gola Dhoro: A reading from animal bones. Antiquity 84, 528–543 (2010).
Google Scholar 

17.
Chase, B. On the pastoral economies of Harappan Gujarat: Faunal analyses at Shikarpur in context. Herit. J. Multidiscip. Stud. Archaeol. 2, 1–22 (2014).
Google Scholar 

18.
Chase, B., Ajithprasad, P., Rajesh, S. V., Patel, A. & Sharma, B. Materializing Harappan identities: Unity and diversity in the borderlands of the Indus Civilization. J. Anthropol. Archaeol. 35, 63–78 (2014).
Google Scholar 

19.
Belcher, W. R. Marine exploitation in the Third Millennium BC- The eastern coast of Pakistan. Paleorient 31, 79–85 (2005).
Google Scholar 

20.
Belcher, W. R. Fish exploitation of the Indus Valley Tradition. In Indus Ethnobiology 95–174 (Lexington Books, Lanham, 2003).
Google Scholar 

21.
Chase, B., Meiggs, D., Ajithprasad, P. & Slater, P. A. What is left behind: Advancing interpretation of pastoral land-use in Harappan Gujarat using herbivore dung to examine bioshphere strontium isotope (87Sr/86Sr) variation. J. Archaeol. Sci. 92, 1–12 (2018).
Google Scholar 

22.
Chase, B., Meiggs, D., Ajithprasad, P. & Slater, P. A. Pastoral land-use of the Indus Civilization in Gujarat: Faunal analyses and biogenic isotopes at Bagasra. J. Archaeol. Sci. 50, 1–15 (2014).
Google Scholar 

23.
Miller, L. J. Secondary products and urbanism in South Asia: The evidence for traction at Harappa. In Indus Ethnobiology: New perspective from the field (eds Weber, S. A. & Belcher, W. R.) 251–326 (Lexington Books, Lanham, 2003).
Google Scholar 

24.
Bourgeois, G. & Gouin, P. Résultats D ’ une analyse de traces organiques fossiles dans une ‘Faisselle’ Harappéenne. Paléorient 21, 125–128 (1995).
Google Scholar 

25.
Evershed, R. P., Dudd, S. N., Copley, M. S. & Mukherjee, A. Identification of animal fats via compound specific δ13C values of individual fatty acids: Assessments of results for reference fats and lipid extracts of archaeological pottery vessels. Doc. Praehist. 29, 73–96 (2002).
Google Scholar 

26.
Copley, M. S., Berstan, R., Straker, V., Payne, S. & Evershed, R. P. Dairying in antiquity. II. Evidence from absorbed lipid residues dating to the British Bronze Age. J. Archaeol. Sci. 32, 505–521 (2005).
Google Scholar 

27.
Spangenberg, J. E., Jacomet, S. & Schibler, J. Chemical analyses of organic residues in archaeological pottery from Arbon Bleiche 3, Switzerland—evidence for dairying in the late Neolithic. J. Archaeol. Sci. 33, 1–13 (2006).
Google Scholar 

28.
Craig, O. E. et al. Stable isotope analysis of Late Upper Palaeolithic human and faunal remains from Grotta del Romito (Cosenza) Italy. J. Archaeol. Sci. 37, 2504–2512 (2010).
Google Scholar 

29.
Craig, O. E. et al. Earliest evidence for the use of pottery. Nature 496, 351–354 (2013).
CAS  PubMed  ADS  Google Scholar 

30.
Craig, O. E., Taylor, G., Mulville, J., Collins, M. J. & Parker, P. M. The identification of prehistoric dairying activities in the Western Isles of Scotland: an integrated biomolecular approach. J. Archaeol. Sci. 32, 91–103 (2005).
Google Scholar 

31.
Evershed, R. P. Organic residue analysis in archaeology: The archaeological biomarker revolution. Archaeometry 50, 895–924 (2008).
CAS  Google Scholar 

32.
Craig, O. E., Love, G. D., Isaksson, S., Taylor, G. & Snape, C. E. Stable carbon isotopic characterisation of free and bound lipid constituents of archaeological ceramic vessels released by solvent extraction, alkaline hydrolysis and catalytic hydropyrolysis. J. Anal. Appl. Pyrolysis 71, 613–634 (2004).
CAS  Google Scholar 

33.
Dudd, S. N. & Evershed, R. P. Direct demonstration of milk as an element of archaeological economies. Science 282, 1478–1481 (1998).
CAS  PubMed  Google Scholar 

34.
Buonasera, T. Y., Tremayne, A. H., Darwent, C. M., Eerkens, J. W. & Mason, O. K. Lipid biomarkers and compound specific δ13C analysis indicate early development of a dual-economic system for the Arctic Small Tool tradition in northern Alaska. J. Archaeol. Sci. 61, 129–138 (2015).
CAS  Google Scholar 

35.
Meier-Augenstein, W. Stable isotope analysis of fatty acids by gas chromatography–isotope ratio mass spectrometry. Anal. Chim. Acta 465, 63–79 (2002).
CAS  Google Scholar 

36.
Mottram, H. R., Dudd, S. N., Lawrence, G. J., Stott, A. W. & Evershed, R. P. New chromatographic, mass spectrometric and stable isotope approaches to the classification of degraded animal fats preserved in archaeological pottery. J. Chromatogr. A 833, 209–221 (1999).
CAS  Google Scholar 

37.
Gregg, M. W., Banning, E. B., Gibbs, K. & Slater, G. F. Subsistence practices and pottery use in Neolithic Jordan: Molecular and isotopic evidence. J. Archaeol. Sci. 36, 937–946 (2009).
Google Scholar 

38.
Copley, M. S. et al. Direct chemical evidence for widespread dairying in prehistoric Britain. Proc. Natl. Acad. Sci. 100, 1524–1529 (2003).
CAS  PubMed  ADS  Google Scholar 

39.
Dunne, J., di Lernia, S., Chłodnicki, M., Kherbouche, F. & Evershed, R. P. Timing and pace of dairying inception and animal husbandry practices across Holocene North Africa. Quat. Int. https://doi.org/10.1016/j.quaint.2017.06.062 (2017).
Article  Google Scholar 

40.
Dunne, J. et al. First dairying in green Saharan Africa in the fifth millennium BC. Nature 486, 390–394 (2012).
CAS  PubMed  ADS  Google Scholar 

41.
Copley, M. S., Clark, K. & Evershed, R. P. Organic-residue analysis of pottery vessels and clay balls. In Changing Materialities at Çatalhoyuk: Reports from the 1995–99 Seasons 169–174 (McDonald Institute for Archaeological Research, Cambridge, 2005).
Google Scholar 

42.
Roffet-Salque, M., Lee, M. R. F., Timpson, A. & Evershed, R. P. Impact of modern cattle feeding practices on milk fatty acid stable carbon isotope compositions emphasise the need for caution in selecting reference animal tissues and products for archaeological investigations. Archaeol. Anthropol. Sci. 9, 1343–1348 (2017).
Google Scholar 

43.
Craig, O. E. et al. Ancient lipids reveal continuity in culinary practices across the transition to agriculture in Northern Europe. Proc. Natl. Acad. Sci. 108, 17910–17915 (2011).
CAS  PubMed  ADS  Google Scholar 

44.
Salque, M. et al. Earliest evidence for cheese making in the sixth millennium BC in northern Europe. Nature 493, 522–525 (2013).
CAS  PubMed  ADS  Google Scholar 

45.
Correa-Ascencio, M., Robertson, I. G., Cabrera-Cortés, O., Cabrera-Castro, R. & Evershed, R. P. Pulque production from fermented agave sap as a dietary supplement in Prehispanic Mesoamerica. Proc. Natl. Acad. Sci. 111, 14223–14228 (2014).
CAS  PubMed  ADS  Google Scholar 

46.
Kimpe, K., Jacobs, P. A. & Waelkens, M. Analysis of oil used in late Roman oil lamps with different mass spectrometric techniques revealed the presence of predominantly olive oil together with traces of animal fat. J. Chromatogr. A 937, 87–95 (2001).
CAS  PubMed  Google Scholar 

47.
Eerkens, J. The preservation and identification of Piñon resins by GC-MS in pottery from the western Great Basin. Archaeometry 44, 95–105 (2002).
CAS  Google Scholar 

48.
Gregg, M. W., Brettell, R. & Stern, B. Bitumen in Neolithic Iran: Biomolecular and isotopic evidence. In Archaeological Chemistry 137–151 (American Chemical Society, Washington, 2007).
Google Scholar 

49.
Lucquin, A., March, R. J. & Cassen, S. Analysis of adhering organic residues of two “coupes-à-socles” from the Neolithic funerary site “La Hougue Bie” in Jersey: evidences of birch bark tar utilisation. J. Archaeol. Sci. 34, 704–710 (2007).
Google Scholar 

50.
Stacey, R., Cartwright, C., Tanimoto, S. & Villing, A. Coatings and contents: Investigations of residues on four fragmentary sixth-century B.C. vessels from Naukratis (Egypt). Br. Museum Tech. Res. Bull. 4, 19–26 (2010).
Google Scholar 

51.
Brecoulaki, H., Andreotti, A., Bonaduce, I., Colombini, M. P. & Lluveras, A. Characterization of organic media in the wall-paintings of the “Palace of Nestor” at Pylos, Greece: Evidence for a secco painting techniques in the Bronze Age. J. Archaeol. Sci. 39, 2866–2876 (2012).
CAS  Google Scholar 

52.
Spades, S. & Russ, J. GC–MS analysis of lipids in prehistoric rock paints and associated oxalate coatings from the Lower Pecos Region, Texas. Archaeometry 47, 115–126 (2005).
CAS  Google Scholar 

53.
Eckmeier, E. & Wiesenberg, G. L. B. Short-chain n-alkanes (C16–20) in ancient soil are useful molecular markers for prehistoric biomass burning. J. Archaeol. Sci. 36, 1590–1596 (2009).
Google Scholar 

54.
Chakraborty, K. S. et al. Enamel isotopic data from the domesticated animals at Kotada Bhadli, Gujarat, reveals specialized animal husbandry during the Indus Civilization. J. Archaeol. Sci. Rep. 21, 2 (2018).
Google Scholar 

55.
Chakraborty, K. S. Subsistence-Based Economy and the Regional Interaction Processes of the Indus Civilization Borderland in Kachchh, Gujarat: A Bio-Molecular Perspective (University of Toronto, Toronto, 2019).
Google Scholar 

56.
Shirvalkar, P. & Rawat, Y. S. Excavation at Kotada Bhadli, District Kachchh, Gujarat: A perliminary report. Puratattva 42, 182–201 (2012).
Google Scholar 

57.
Goyal, P. Observations on faunal remains recovered from Kotada Bhadli. In Excavation at Kotada Bhadli (eds Shirvalkar, P. & Prasad, E.) 136–149 (Archaeological Survey of India, New Delhi, 2020).
Google Scholar 

58.
Joglekar, P. P. & Goyal, P. Faunal remains from Shikarpur, a Harappan site in Gujarat India. Iran. J. Archaeol. Stud. 1, 15–25 (2011).
Google Scholar 

59.
Goyal, P. & Joglekar, P. P. Archaeozoological remains from the site of Kanmer. In Excavation at Kanmer (2005–2006 to 2008–2009): Kanmer archaeological research project an Indo-Japanese collaboration (eds Kharakwal, J. S. et al.) 767–794 (Indus Project Research Institute for Humanity and Nature, Kyoto, 2012).
Google Scholar 

60.
Meadow, R. H. & Patel, A. K. Prehistoric pastoralism in northwestern South Asia from the Neolithic through the Harappan Period. In Indus Ethnobiology: New perspective from the field (eds Weber, S. & Belcher, W. R.) 65–94 (Lexington Books, Lanham, 2003).
Google Scholar 

61.
Patel, A. The Primary Pastoral economy of Dholavira: A first look at animals and urban life in third millennium Kutch. In South Asian Archaeology 1995 (ed. Allchin, B.) 101–113 (The Ancient India and Iran Trust, Cambridge, 1997).
Google Scholar 

62.
Miller, L. J. Urban Economies in Early States: The Secondary Products Revolution in the Indus Civilization (New York University, New York, 2004).
Google Scholar 

63.
Halstead, P. Mortality models and milking: Problems of uniformitarianism, optimality and equifinality reconsidered. Anthropozoologica https://doi.org/10.4319/lo.2013.58.2.0489 (1998).
Article  Google Scholar 

64.
Gillis, R. E. et al. The evolution of dual meat and milk cattle husbandry in Linearbandkeramik societies. Proc. R. Soc. B Biol. Sci. 284, 2 (2017).
Google Scholar 

65.
Sternberg, L. O., Deniro, M. J. & Johnson, H. B. Isotope ratios of cellulose from plants having different photosynthetic pathways. Plant Physiol. 74, 557–561 (1984).
CAS  PubMed  PubMed Central  Google Scholar 

66.
Zhang, C. et al. Diets and environments of late Cenozoic mammals in the Qaidam Basin, Tibetan Plateau: Evidence from stable isotopes. Earth Planet. Sci. Lett. 333–334, 70–82 (2012).
ADS  Google Scholar 

67.
Cerling, T. E. et al. Global vegetation change through the Miocene/Pliocene boundary. Nature 389, 153–158 (1997).
CAS  ADS  Google Scholar 

68.
Sternberg, L. O., Deniro, M. J. & Ting, I. P. Carbon, hydrogen, and oxygen isotope ratios of cellulose from plants having intermediary photosynthetic modes. Plant Physiol. 74, 104–107 (1984).
CAS  PubMed  PubMed Central  Google Scholar 

69.
Smith, B. N. & Epstein, S. Two categories of 13C/12C ratios for higher plants. Plant Physiol. 47, 380–384 (1971).
CAS  PubMed  PubMed Central  Google Scholar 

70.
Correa-Ascencio, M. & Evershed, P. R. High throughput screening of organic residues in archaeological potsherds using direct acidified methanol extraction. Anal. Methods 6, 1330–1340 (2014).
CAS  Google Scholar 

71.
Evershed, R. P., Heron, C. & Goad, L. J. Analysis of organic residues of archaeological origin by high-temperature gas chromatography and gas chromatography-mass spectrometry. Analyst 115, 1339–1342 (1990).
CAS  ADS  Google Scholar 

72.
Papakosta, V., Smittenberg, R. H., Gibbs, K., Jordan, P. & Isaksson, S. Extraction and derivatization of absorbed lipid residues from very small and very old samples of ceramic potsherds for molecular analysis by gas chromatography-mass spectrometry (GC-MS) and single compound stable carbon isotope analysis by gas chromatogra. Microchem. J. 123, 196–200 (2015).
CAS  Google Scholar 

73.
Demirci, Ö, Lucquin, A., Craig, O. E. & Raemaekers, D. C. M. First lipid residue analysis of Early Neolithic pottery from Swifterbant ( the Netherlands, ca 4300–4000 BC ). Archaeol. Anthropol. Sci. https://doi.org/10.1007/s12520-020-01062-w (2020).
Article  Google Scholar 

74.
Carrer, F. et al. Chemical analysis of pottery demonstrates prehistoric origin for high-altitude alpine dairying. PLoS ONE 11, e0151442 (2016).
PubMed  PubMed Central  Google Scholar 

75.
Heron, C. et al. First molecular and isotopic evidence of millet processing in prehistoric pottery vessels. Sci. Rep. 6, 1–9 (2016).
Google Scholar 

76.
Pecci, A. & Cau Ontiveros, M. A. Report on the Analyses of the Organic Residues in Archaeological Samples from the Project ‘Excavating the Roman Peasant’. University of Barcelona (2010).

77.
Gregg, M. W. & Slater, G. F. A new method for extraction, isolation and transesterification of free fatty acids from archaeological pottery. Archaeometry 52, 833–854 (2010).
CAS  Google Scholar 

78.
Copley, M. S. et al. Detection of palm fruit lipids in archaeological pottery from Qasr Ibrim, Egyptian Nubia. Proc. R. Soc. B Biol. Sci. 268, 593–597 (2001).
CAS  Google Scholar 

79.
Evershed, R. P., Copley, M. S., Dickson, L. & Hansel, F. A. Experimental evidence for the processing of marine animal products and other commodities containing polyunsaturated fatty acids in pottery vessels. Archaeometry 50, 101–113 (2008).
CAS  Google Scholar 

80.
Hansel, F. A., Copley, M. S., Madureira, L. A. S. & Evershed, R. P. Thermally produced ω-(o-alkylphenyl)alkanoic acids provide evidence for the processing of marine products in archaeological pottery vessels. Tetrahedron Lett. 45, 2999–3002 (2004).
CAS  Google Scholar 

81.
Meadow, R. H. Prehistoric wild sheep and sheep domestication on the eastern margin of the Middle East. in Animal Domestication and Its Cultural Context (eds. Crabtree, P. J., Campana, D. V. & Ryan, K.) 24–36 (University Museum, Univerity of Pennsylvania, 1989).

82.
Craig, O. E. et al. Distinguishing wild ruminant lipids by gas chromatography/combustion/isotope ratio mass spectrometry. Rapid Commun. Mass Spectrom. 26, 2359–2364 (2012).
CAS  PubMed  ADS  Google Scholar 

83.
Pokharia, A. K. Floral Remains. in Excavation at Kanmer (2005–06 to 2008–09): Kanmer archaeological research project an Indo-Japanese collaboration (eds. Kharakwal, J. S., Rawat, Y. S. & Osada, T.) 795–812 (Indus project, research institute for humanity and Nature, 2012).

84.
Steele, V. J., Stern, B. & Stott, A. W. Olive oil or lard?: distinguishing plant oils from animal fats in the archeological record of the eastern Mediterranean using gas chromatography/combustion/isotope ratio mass spectrometry. Rapid Commun. Mass Spectrom. 24, 3478–3484 (2010).
CAS  PubMed  ADS  Google Scholar 

85.
Spangenberg, J. E. & Ogrinc, N. Authentication of vegetable oils by bulk and molecular carbon isotope analyses with emphasis on olive oil and pumpkin seed oil. J. Agric. Food Chem. 49, 1534–1540 (2001).
CAS  PubMed  Google Scholar 

86.
Hammann, S. & Cramp, L. J. E. Towards the detection of dietary cereal processing through absorbed lipid biomarkers in archaeological pottery. J. Archaeol. Sci. 93, 74–81 (2018).
CAS  Google Scholar 

87.
Colonese, A. C. et al. New criteria for the molecular identification of cereal grains associated with archaeological artefacts. Sci. Rep. 7, 1–7 (2017).
CAS  Google Scholar 

88.
Courel, B. et al. Organic residue analysis shows sub-regional patterns in the use of pottery by Northern European hunter-gatherers. R. Soc. Open Sci. 7, 2 (2020).
Google Scholar 

89.
Hendy, J. et al. Ancient proteins from ceramic vessels at Çatalhöyük West reveal the hidden cuisine of early farmers. Nat. Commun. 9, 2 (2018).
Google Scholar 

90.
Chase, B., Meiggs, D. & Ajithprasad, P. Pastoralism, climate change, and the transformation of the Indus Civilization in Gujarat: Faunal analyses and biogenic isotopes. J. Anthropol. Archaeol. 59, 101173 (2020).
Google Scholar 

91.
Margabandhu, C. Technology of trasport vehicles in early India. In Radiocarbon and Indian Archaeology (eds Agarwal, D. P. & Ghosh, A.) 182–189 (Munshiram Manoharlal Publishers, New Delhi, 1973).
Google Scholar 

92.
Fairservis, W. J. Cattle. Exped. Mag. 28, 43–50 (1986).
Google Scholar 

93.
Chase, B. Family matters in Gujarat. In Walking with the Unicorn: Social Organization and Material Culture in Ancient South Asia (eds Frenez, D. et al.) 90–103 (Archaeopress Publishing LTD, Oxford, 2018).
Google Scholar  More

1.
Makino, S. Biology of Latibulus argiolus (Hymenoptera, Ichneumonidae), a Parasitoid of the Paper Wasp Polistes biglumis (Hymenoptera, Vespidae). Kontyû 51, 426–434 (1983).
Google Scholar 
2.
Oh, S.-H., An, S.-L. & Lee, J.-W. Review of Korean Latibulus (Hymenoptera: Ichneumonidae: Cryptinae) and a key to the world species. Can. Entomol. 144, 509–525 (2012).
Article  Google Scholar 

3.
Quicke, D. L. J. The Braconid and Ichneumonid Parasitoid Wasps: Biology Systematics, Evolution and Ecology (Wiley Blackwell, Amsterdam, 2015).
Google Scholar 

4.
Paukku, S. & Kotiaho, J. S. Female oviposition decisions and their impact on progeny life-history traits. J. Insect Behav. 21, 505–520 (2008).
Article  Google Scholar 

5.
Rusina, L. Y. The role of parasitoids in regulation of Polistes Wasp population (Hymenoptera, Vespidae: Polistinae). Entomol. Rev. 93, 271–280 (2012).
Article  Google Scholar 

6.
Coelho, N. H. P. et al. Understanding genetic diversity, spatial genetic structure, and mating system through microsatellite markers for the conservation and sustainable use of Acrocomia aculeata (Jacq.) Lodd. Ex Mart. Conserv. Genet. 19, 879–891 (2018).
CAS  Article  Google Scholar 

7.
Manlik, O. et al. Demography and genetics suggest reversal of dolphin source-sink dynamics, with implications for conservation. Mar. Mammal Sci. 35, 732–759 (2019).
Article  Google Scholar 

8.
Nowicki, P. et al. What keeps “living dead” alive: demography of a small and isolated population of Maculinea (=Phengaris) alcon. J. Insect Conserv. 23, 201–210 (2019).
Article  Google Scholar 

9.
Selkoe, K. A. & Toonen, R. J. Microsatellites for ecologists: a practical guide to using and evaluating microsatellite markers. Ecol. Lett. 9, 615–629 (2006).
Article  Google Scholar 

10.
Olafsson, K., Hjorleifsdottir, S., Pampoulie, C., Hreggvidsson, G. O. & Gudjonsson, S. Novel set of multiplex assays (SalPrint15) for efficient analysis of 15 microsatellite loci of contemporary samples of the Atlantic salmon (Salmo salar). Mol. Ecol. Resour. 10, 533–537 (2010).
CAS  Article  Google Scholar 

11.
Randi, E. et al. Multilocus detection of wolf x dog hybridization in Italy, and guidelines for marker selection. PLoS ONE 9, e86409 (2014).
ADS  Article  Google Scholar 

12.
Hung, C.-M., Yu, A.-Y., Lai, Y.-T. & Shaner, P.-J.L. Developing informative microsatellite markers for nonmodel species using reference mapping against a model species’ genome. Sci. Rep. 6, 23087 (2016).
ADS  CAS  Article  Google Scholar 

13.
Marcus, T., Assmann, T., Durka, W. & Drees, C. A suite of multiplexed microsatellite loci for the ground beetle Abax parallelepipedus (Piller and Mitterpacher, 1783) (Coleoptera, Carabidae). Conserv. Genet. Resour. 5, 1151–1156 (2013).
Article  Google Scholar 

14.
Panagiotopoulou, H., Baca, M., Baca, K., Stanković, A. & Żmihorski, M. Optimization and validation of a multiplex assay for microsatellite loci analysis in the field cricket, Gryllus campestris (Orthoptera: Gryllidae). J. Asia-Pac. Entomol. 18, 421–424 (2015).
CAS  Article  Google Scholar 

15.
PacBio Pacific Biosciences, Procedure & Checklist: 2 kb Template Preparation and Sequencing. https://www.pacb.com/wp-content/uploads/2015/09/Procedure-Checklist-2-kb-Template-Preparation-and-Sequencing.pdf.

16.
Faircloth, B. C. MSATCOMMANDER: detection of microsatellite repeat arrays and automated, locus-specific primer design. Mol. Ecol. Resour. 8, 92–94 (2008).
CAS  Article  Google Scholar 

17.
Zhang, Z., Schwartz, S., Wagner, L. & Webb, M. A greedy algorithm for aligning DNA sequences. J. Comput. Biol. 7, 203–214 (2000).
CAS  Article  Google Scholar 

18.
Morgulis, A. et al. Database indexing for production MegaBLAST searches. Bioinformatics 24, 1757–1764 (2008).
CAS  Article  Google Scholar 

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

20.
Schuelke, M. An economic method for the fluorescent labeling of PCR fragments. Nat. Biotechnol. 18, 233–234 (2000).
CAS  Article  Google Scholar 

21.
Austin, J. D. et al. Permanent genetic resources added to Molecular Ecology Resources Database 1 February 2011–31 March 2011. Mol. Ecol. Resour. 11, 757–758 (2011).
CAS  Article  Google Scholar 

22.
Peakall, R. & Smouse, P. E. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics 28, 2537–2539 (2012).
CAS  Article  Google Scholar 

23.
Van Oosterhout, C., Hutchinson, W. F., Wills, D. P. M. & Shipley, P. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 4, 535–538 (2004).
Article  Google Scholar 

24.
Kalinowski, S. T., Taper, M. L. & Marshall, T. C. Revising how the computer program CERVUS accommodates genotyping error increases success in paternity assignment. Mol. Ecol. 16, 1099–1106 (2007).
Article  Google Scholar 

25.
Rousset, F. GENEPOP’007: a complete reimplementation of the GENEPOP software for Windows and Linux. Mol. Ecol. Resour. 8, 103–106 (2008).
Article  Google Scholar 

26.
Rice, W. R. Analyzing tables of statistical tests. Evolution 43, 223–225 (1989).
Article  Google Scholar 

27.
Pompanon, F., Bonin, A., Bellemain, E. & Taberlet, P. Genotyping errors: causes, consequences and solutions. Nat. Rev. Genet. 6, 846–847 (2005).
Article  Google Scholar 

28.
Grohme, M. A., Soler, R. F., Wink, M. & Frohme, M. Microsatellite marker discovery using single molecule real-time circular consensus sequencing on the Pacific Biosciences RS. Biotechniques 55, 253–256 (2013).
CAS  Article  Google Scholar 

29.
Liljegren, M. M., de Muinck, E. J. & Trosvik, P. Microsatellite length scoring by single molecule real time sequencing-effects of sequence structure and PCR regime. PLoS ONE 11, e0159232 (2016).
Article  Google Scholar 

30.
Dutta, N. et al. Microsatellite marker set for genetic diversity assessment of primitive Chitala chitala (Hamilton, 1822) derived through SMRT sequencing technology. Mol. Biol. Rep. 46, 41–49 (2018).
Article  Google Scholar 

31.
Wei, N., Bemmels, J. B. & Dick, C. W. The effects of read length, quality and quantity on microsatellite discovery and primer development: from Illumina to Pac Bio. Mol. Ecol. Resour. 14, 953–965 (2014).
CAS  PubMed  Google Scholar 

32.
Corner, S., Yuzbasiyan-Gurkan, V., Agnew, D. & Venta, P. J. Development of a 12-plex of new microsatellite markers using a novel universal primer method to evaluate the genetic diversity of jaguars (Panthera onca) from North American zoological institutions. Conservation Genet. Resour. 11, 487–497 (2019).
Article  Google Scholar 

33.
Chapuis, M.-P. & Estoup, A. Microsatellite null alleles and estimation of population differentiation. Mol. Biol. Evol. 24, 621–631 (2007).
CAS  Article  Google Scholar 

34.
Blondin, L. et al. Characterization and comparison of microsatellite markers derived from genomic and expressed libraries for the desert locust. J. Appl. Entomol. 137, 673–683 (2013).
Article  Google Scholar 

35.
Chistiakov, D. A., Hellemans, B. & Volckaert, F. A. Microsatellites and their genomic distribution, evolution, function and applications: a review with special reference to fish genetics. Aquaculture 255, 1–29 (2006).
CAS  Article  Google Scholar 

36.
Cheng, L., Zhang, Y., Lu, C.-Y., Li, C. & Sun, X.-W. Development and characterization of four moderate multiplex microsatellite panels in crucian carp (Carassius auratus). Conserv. Genet. Resour. 5, 821–823 (2013).
Article  Google Scholar  More

To investigate the pathway from adult population characteristics to spawning behaviour, egg production, and ultimately to recruitment (Fig. 1), we used three data sources; an egg survey (for estimates of egg distribution, total egg production, and environmental variables), biological samples of the commercial fishery (for estimates of spawning duration and peak, and maternal body condition), and stock assessment outputs (for estimates of age-1 recruits, spawning stock biomass and age structure).
Figure 1

Conceptual framework of the pathway from spawners to recruits and the underlying mechanisms investigated (stock demographic structure and environmental conditions in red and green, respectively).

Full size image

Egg survey data
Sampling
Mackerel enter the southern Gulf of St. Lawrence (sGSL, Eastern Canada) in early June each year to spawn, after overwintering along the north-eastern US continental shelf (from Sable Island to the Mid-Atlantic Bight29,30). Each year, Fisheries and Oceans Canada (DFO) conducts a 2-week long mackerel egg survey in the sGSL (a 65-station fixed grid 20 nautical miles apart spanning the dominant mackerel spawning area) around the average mackerel peak spawning date of June 21st. Over this period, a large fraction of spawning occurs and the survey is therefore believed to reflect appropriately spawning intensity and spatio-temporal properties. Stations consist of double oblique tows using 61-cm Bongo nets with 333 µm mesh size and flowmeters carried out on board a research vessel at a speed of 2.5 knots from 0 to 50 m depth to estimate daily and total egg production while also measuring physical and biological oceanographic variables (see further details in SI Appendix A). This survey has been carried out consistently since 1982, except for no surveys in 1995 and 1997. Several indices are derived from this mackerel egg survey: total egg production, egg distribution, water temperature, and zooplankton biomass, species composition, abundance, and distribution.
Total egg production and distribution
Annual total egg production was calculated according to a standard DFO protocol based on the Daily Egg Production Method31. Stage 1 (spawned less than 24 h ago) and 5 (i.e., damaged stage 1 eggs) egg counts were standardized by the volume of filtered water and the depth of the sampled water column to provide egg densities per station (number m−2). These numbers were then adjusted for incubation time32 to obtain daily egg production point estimates. Spatial interpolation was done across a grid of 3320 coordinates using ordinary kriging to calculate a mean daily egg production estimate per grid cell, which was extrapolated to the surface area sampled. Annual egg production estimates were obtained by dividing by the proportion of reproductively active fish at the median date of the survey. This latter value, along with peak spawning date and spawning duration was calculated using a logistic model describing the daily evolution of the gonadosomatic index, based on corresponding biological data (see further details in Doniol-Valcroze et al.31, and in “Commercial fishery sampling”).
To examine the potential inter-annual spatial mismatch between spawning location and the optimal habitat for larvae, we calculated the spatial extent (spawning area) and the position of the centre of gravity (spawning longitude and latitude) of spawning for each year in the time series. The spatial extent of egg production was determined using an α-convex hull on stations where eggs were present33. The centre of gravity of total egg production was calculated by taking the arithmetic mean of the coordinates of each station weighted by their individual observed egg production.
Environmental indices
Sea surface temperature (SST, °C) directly affects early life stage growth and survival7, but might also have an indirect effect on recruitment through adult spawning behaviour, as mackerel generally spawn between 8 and 15 °C34. Therefore, we produced an SST index by averaging June CTD-measured mean water temperatures in the first 10 m over stations, where the majority of mackerel eggs and larvae occur35.
We hypothesized that the main adult mackerel prey (i.e., C. hyperboreus and capelin, Mallotus villosus36) might be influential as well, as they may affect spawning location and therefore be an indirect driver of recruitment. Capelin is despite its importance as prey in terms of weight36 not considered as a potential driver of spawning location, because its consumption by mackerel is infrequent, only important to the larger mackerel and likely opportunistic. As such, habitat selection is most likely to be related to copepod abundance and we developed spatial, biomass, and composition indices in June in the sGSL only for C. hyperboreus. As a proxy of adult mackerel prey location, we computed the annual centre of gravity of C. hyperboreus biomass (latitude and longitude) with the same methodology used for total egg production. Also, we estimated the total C. hyperboreus biomass (mg m−2) in the sGSL37. The percentage of C. hyperboreus biomass relative to the total Calanus spp. biomass (% C. hyp.) was calculated as we hypothesized that changes in C. hyperboreus proportion may have influenced adult mackerel feeding behaviour and thus spawning locations.
Mackerel larvae mainly feed on the early life stages (eggs, nauplii, and young copepodites) of C. finmarchicus, Pseudocalanus spp. and Temora longicornis25. The copepod daily egg production (CEDP, µg egg carbon L−1 d−1) of these three copepod taxa, calculated based on adult female abundance and species-specific per capita daily egg production (see details in the SI Appendix A), was previously recognized as a good predictor of mackerel recruitment23,24,25. High larval prey abundance might, however, be irrelevant when there is a temporal or spatial mismatch with larval distribution. An annual (y) index of a temporal match was therefore calculated in June in the spawning area as the proportion of older stage 6 female C. finmarchicus, producing prey for mackerel early life stages, with respect to the number of younger immature copepodite stages 4 and 526 (Eq. 1).

$${Temporal match}_{y}=100%times {N}_{C. fin female}/{N}_{C. fin stages 4-5}$$
(1)

Higher percentages of stage 6 female copepodites during mackerel spawning (i.e., a later development of the plankton community) should improve the temporal match between hatching and the availability of prey for emerging larvae26. This same index could not include Pseudocalanus spp. and Temora longicornis as only data for stage 6 adults were available. C. finmarchicus is, however, considered to be a good indicator of the overall zooplankton phenology in spring and early summer in the sGSL and should also reflect Pseudocalanus spp. and Temora longicornis phenology27. An annual index of a spatial match between mackerel egg distribution and their near-future prey was determined as the sum of mackerel daily egg production (DEP) at stations (s) with sufficient prey (i.e., copepod daily egg production above a threshold value) divided by the daily egg production of mackerel over all stations (Eq. 2).

$${Spatial match}_{y}=100%times {sum }_{s=1}^{S >threshold}{DEP}_{s,y}/{sum }_{s=1}^{S}{DEP}_{s,y}$$
(2)

The threshold copepod daily egg production value was determined as the 25th quantile of values measured for all years and stations, which excludes zero and near-zero prey availabilities unlikely to be able to support larval survival. This index of spatial match captures a combined effect of the abundance and distribution of the prey in relation to the distribution of the fish eggs. Note that due to the availability of taxonomic zooplankton data, Pseudocalanus spp., Temora spp., C. finmarchicus and C. hyperboreus data and hence all indices derived from it were available for only 21 years (but covering the entire span of the time series; 1982, 1985, 1987, 1990, 1993, 1996, 1999, 2000, 2003 and 2006 to 2017). Spatial and temporal match–mismatch proxies were based on a match with the mackerel eggs rather than the early larval phase. We expect this to introduce little noise as the development time of mackerel eggs is typically less than 6 days and mackerel larval development is fast (about 20 days32). All the environmental variables used and the associated hypotheses are summarized in Table 1.
Table 1 Summary of all the hypotheses tested along the pathway from spawners to recruits and associated references.
Full size table

Commercial fishery sampling
Adult mackerel samples are collected annually by DFO from the commercial fishery. The sampling covers the entire spawning area and period (thrice a week) and on average 4998 (range 421–14,858) individual fish are analysed each year. We used this data to calculate the annual peak spawning date (spawn. peak), spawning duration (spawn. duration), and maternal body condition.
Peak spawning date and duration were calculated each year based on the fit of a logistic model of the daily evolution of the gonadosomatic index. The mean value of the derived symmetrical probability density function was defined as the peak spawning day and the time between the 2.5% and 97.5% quantiles was estimated to represent the spawning duration in days.
As relatively fatter individuals might spawn more and higher quality eggs38, mature females (i.e., reproductive stages 3–839) sampled between their arrival in the sGSL and June 21st (the average peak spawning date) were selected to investigate the potential influence of pre-spawning fat reserves on total egg production and recruitment with the relative body condition index (Kn40, Eq. 3):

$${K}_{n}=frac{W}{{W}_{r}}$$
(3)

where W is the observed somatic weight (g) of an individual and Wr the predicted weight of an individual of a given fork length (FL, cm) calculated with Wr = αFLβ (α and β are nonlinear least-squares regression parameters).
Mackerel SSB, recruitment and age structure
Annual mackerel SSB, recruitment residuals and an index of age structure were derived from an age-structured state-space stock assessment model applied to the period 1968–201828. Note that the model was calibrated using an SSB index directly calculated from total egg production. In the assessment model, a two-parameter Beverton-Holt stock-recruitment relationship was used to estimate annual recruitment (abundance at age 1), and the residuals of this relationship were used in subsequent analyses (Rres). An indicator of the annual age structure was considered as bigger, older mackerel spawners ( > age 5) are known to have a greater fecundity, and spawn in different spatial and temporal niches than younger females35,41. Mean biomass-weighted age (MA) was calculated using mature biomass-at-age (({SSB}_{a})) as follow in the Eq. (4):

$$MA=frac{sum_{a=1}^{a=10}(a{SSB}_{a})}{sum_{a=1}^{a=10}{SSB}_{a}}$$
(4)

MA was based on biomass rather than abundance to better reflect the stock’s reproductive potential42.
Mackerel early life stages are prey for pelagic fish sharing the surface waters of the sGSL. Herring are, relative to other potential predators, dominant, widely distributed and known predators of mackerel eggs and larvae36. Hence, we used cumulated spring and fall herring model-derived annual biomass43 as a proxy of predation pressure on mackerel early life stages.
Statistical analyses
Recruitment variability driven by spawning aspects and environmental gradients
We analysed the relationships between the successive steps leading to recruitment (spawning aspects, egg production and recruitment) and both demographic and environmental effects using generalised linear models (GLMs). All model configurations (response and explanatory variables) are given in Supplementary Table S2. Explanatory variables were normalized (i.e., by subtracting the mean and dividing by the standard deviation for each variable) to facilitate comparison of their respective effects (i.e., through their coefficients). When the response variable was Rres (with a 1-year lag), residuals were assumed to follow a Gaussian distribution with an identity link function, whereas for the other response variables a Gamma distribution with a log link function was used (as they can only take positive values44). Before performing GLM computations, collinearity between explanatory variables was measured using variance inflation factors (VIFs), considering a VIF threshold of 344. Specifically, mackerel SSB and MA were highly correlated (Pearson correlation coefficient  > 0.7, see Supplementary Fig S1), so distinct sets of GLMs testing SSB or MA on spawning aspects were used. A backwards model selection procedure was performed, choosing the model with the lowest Akaike’s information criterion corrected for small samples sizes (AICc). If independent models including either SSB or MA showed an AICc difference less than 2, both were reported. Assumptions of homoscedasticity and normality were checked using residual plots while assumptions of independence (to ensure no autocorrelation was present) were checked using correlograms. By replacing GLMs with generalized additive models, the same conclusions were reached and there were no indications of strong non-linear effects.
Variability in total egg production (TEP) could not be linked directly to SSB and MA using regression techniques, because of model circularity (a TEP derived SSB index was used to estimate SSB) and collinearity (SSB and MA are significantly correlated and difficult to disentangle). Although the relative effect size of both variables could not be measured, the positive link between them is well established in the literature (i.e., that larger, older fish produce more eggs41). We, therefore, focussed our efforts on the possible link between TEP per unit of biomass, thereby removing the effect of fish number- and weight-at-age, and maternal body condition. Furthermore, by working with stock–recruitment residuals, we removed in large part the intrinsically related process of TEP. That is, the stock–recruitment relationship is presumably created by the biological dependence of TEP on SSB, and subsequently of recruitment on TEP. This link was hence not explicitly considered, although being present. A Jackknife procedure was conducted to assess the consistency and robustness of the optimal models explaining recruitment residuals (see SI appendix A). Also, recruitment estimates are inherently dependent on the modelling choices45, and we verified that recruitment residuals obtained under different assumptions (i.e., through a Virtual Population Analysis, VPA46) were not differently explained by the considered variables (see SI appendix A for more details).
Stability of the recruitment-larval prey availability relationship
Since Castonguay et al.23, a different stock assessment model has been employed, resulting in new recruitment timeseries47. As a baseline for comparison, we, therefore, refitted the recruitment–CEDP relationship from Castonguay et al.23 with the updated estimates and including all years (1982–2017, linear modelling). We hypothesized that, with the addition of new years of data, potential changes in the performance of this quantitative food index (i.e., CEDP) in predicting recruitment would be driven by a temporal change in the relationship because of altering underlying mechanisms. The latter could manifest itself as changes in the spatial or temporal match between the CEDP and the spawning distribution (a proxy of larval distribution), i.e., the ‘effective’ prey availability. Thus, we examined whether changing larval prey availability in space and time, coupled with a changing mackerel larval quality (using adult Kn as a proxy), can explain residuals and the potential breakdown of the Rres-CEDP relationship. Then, the drivers behind the spatial match-mismatch between mackerel eggs and larval prey were investigated. We considered maternal body condition, SST, and C. hyperboreus longitude (i.e., spawner prey). We also retained the relative abundance of C. hyperboreus in the Calanus spp. community (% C. hyp.), as this species does not produce eggs and nauplii available to mackerel larvae in the summer in the sGSL37,48 and appears to reduce abundance of C. finmarchicus early life stages (i.e., mackerel larval prey) through predation49. Thus, years with a large proportion of C. hyperboreus in the plankton community may display a larger mismatch between mackerel eggs and CEDP. A beta regression model was used to study the spatial match (as it is a proportion). All statistical analyses were conducted with R (version 3.3.250).
Ethical approval
This study was approved by DFO Research Ethics Board and conducted with methods in accordance with the Canadian Council on Animal Care (ISBN: 0-919087-43-4). More

Mobile phone data can be used to inform different aspects of COVID-19 response (Table 1). At the population level, quantifying changes in human mobility or clustering can help evaluate the impact of an NPI and identify hotspots where additional or different interventions may need to be applied. At the individual level, mobile phone data may be used to understand patterns of individual contacts and enhance contact tracing.
Table 1 Summary of types, metrics, and proposed applications of mobile phone data.
Full size table

Evaluating current interventions and monitoring their release
The most widely used application of mobile phone data in public health to date is the use of telecom geolocation data to track population movements11,12. Mobile phone operators routinely collect Call Detail Records (CDRs) that contain a timestamp and GPS location with a unique identifier for all subscribers. These data thus are typically readily available and offer high coverage to estimate mobility patterns of individuals using their mobile devices. We note that similar time-resolved GPS location data may be passively collected through certain applications, though typically for only a subset of subscribers that may introduce further bias.
CDRs can be used to generate a number of metrics for characterizing large, population-level mobility patterns. Origin-Destination (OD) matrices reflect the number of times a trip is made between two locations (of varying spatial resolution) in a certain period. These matrices can be analyzed over time to detect temporal trends (i.e., holidays, seasonality, weekday vs weekend) and regular hotspots of attraction. These spatial and temporal flows of individuals between locations, including the magnitude and frequency of these movements, can be used to understand the risk of importation from areas with ongoing outbreaks to areas without sustained transmission where there is a risk of reintroduction and resurgence. Aggregate flows can also be used to retrace the likely introduction and spread of an outbreak in new areas and to inform future projections of disease risk or burden across space and decision making around the design and implementation of travel restrictions or increased surveillance.
Aggregate mobility patterns may also be critical pieces of evidence when evaluating the effectiveness of various NPIs. Most NPIs are reliant on modifying physical behavior. Monitoring the volume, frequency, and average distance of flow during interventions can be used to directly quantify the adoption and effect of these interventions, and identify areas of high potential risk to target with different interventions. There are already identified associations between reductions in population-level mobility within and between different locations and COVID-19 incidence6,10,29, though further exploration of which population-level metrics are most closely related to changes in disease risk and whether these associations are sustained throughout an outbreak is needed30. These associations would ideally be interrogated to identify individual behaviors associated with mobility measures that are also associated with individual risk of COVID-19.
The effect on NPIs can also be monitored through subscriber density metrics that combine the recorded GPS location and timestamp of CDRs to capture the real-time population density and identify potential hotspots. When using finer-scale GPS location data, these density metrics may quantify the likelihood or frequency that users came into proximal contact. A third metric derived from CDR or GPS location data, the radius of gyration, quantifies the range over which a single person may travel in a specified time period. Importantly, the data required for these applications are non-identifiable; they cannot be used to identify any given individual’s interactions, but provide population-level insight into the average clustering and movement of individuals. These metrics, along with traditional OD matrix flows, were recently employed in Italy as a way to evaluate the impact of its national lockdown31. Traffic flow between provinces and probability of colocation were reduced initially in the northern provinces, where the COVID-19 outbreak was first observed, a clear signal of reactive social distancing. As the epidemic progressed, and especially once the national lockdown was enforced, the entire country saw a reduction in traffic between provinces; however, the probability of colocation remained highly dependent on province and was likely attributed to the number of cases reported in each province. Interestingly, the average distance traveled by individuals was significantly reduced across all provinces after the initial outbreak was confirmed.
The use of Bluetooth data (records of proximal interactions between Bluetooth-enabled devices) to quantify physical clustering or real-time density of subscribers at small spatial scales (e.g., zip codes) and fine temporal resolution has been explored for the purposes of contact tracing (see below). The use of these data has been considered less for population-level analyses, though it offers another source of information on behavioral changes under different NPIs. When activated, mobile phones will emit a Bluetooth beacon that is detected by other activated phones. When two Bluetooth-enabled devices are within range, the date, time, distance and duration of interaction can be recorded. The frequency or number of these interactions (analyzed anonymously to form, broadly, measures of clustering or proximal interaction rates over time) may be important given the role of sustained interaction or overcrowding of individuals32,33,34 and contact structure in SARS-CoV-2 transmission35. Furthermore, Bluetooth data in combination with GPS data or a network of Bluetooth sensors can be used to quantify the amount of time people spend at home or other identified locations when lockdown measures are in place to determine if policies are effective.
These data and measures of population-level mobility or clustering patterns would be exceedingly difficult to collect on a similar scale without mobile phone data. These data are often continuously collected, in near real-time, allowing for continued analysis as an outbreak unfolds. Importantly, though, a baseline understanding of contact or clustering patterns prior to any interventions is necessary to inform estimates of intervention impact.
Facilitating contact tracing
Opt-in applications (apps)36,37,38,39,40,41,42 that rely on digital approaches to enumerate and contact individuals who may have been in proximity with someone infected with COVID-19 have been proposed to increase efficiency and decrease the very large burden of manual contact tracing programs43,44,45. By enabling rapid tracing of perhaps higher proportions of affected individuals, these apps can reduce the amount of time that a potentially infected person would have to infect others, particularly in asymptomatic or pre-symptomatic phases of infection46. Most contact tracing apps collect Bluetooth and/or GPS location data to create trails of contacts over a moving time window (14-28 days). Unlike the data needed to understand population-level, aggregated behaviors described above, these data must be linked to single individuals and capture pairwise interactions with other identifiable individuals. Once a case has been identified, they are added to a list of infected users that is queried by the other phones in the network. If the infected user is detected in the trail of contacts, then the user and their contacts are alerted, either by the app or by a public health official, to initiate isolation and quarantine.
This contact tracing process occurs either in a centralized manner, where user information is sent to a remote computer where matching occurs, or in a decentralized manner, where the matching process occurs on the user’s phone. In order for these approaches to feed directly into public health decision making, a direct line between the developers, public health response teams, and users needs to be put in place. This will also be key to mitigating any privacy concerns, which should be dealt with in a transparent and direct manner. Although there has been little discussion to date, routinely collected, individually-identifiable Bluetooth or fine-scale GPS location data may also be used to infer and quantify high-resolution proximity network structures which may further inform contact tracing efforts, but will also raise additional privacy concerns47,48.
Frameworks to process and analyze mobile phone data
Luckily, computing resources and methods to analyze and extract these data will not likely be the limiting factor in these instances. Groups such as Flowminder and Telenor Research Group have worked for multiple years to develop more streamlined processes to analyze these data, particularly aggregate mobility data, that are able to directly interface with mobile phone operators. Flowminder has produced a suite of CDR aggregates, such as counts of active subscribers per region or counts of travelers, that can then be used to calculate indicators of mobility, such as crowdedness, population mixing, locations of interest, and intra-/inter-regional travel49. The code to extract these metrics is publicly available at50. Telenor Research Group works directly with mobile phone operators to provide researchers with spatially aggregated CDR/mobility data51. Facebook’s Data For Good program provides aggregated mobility data to researchers that come from their subscribers, and companies like Cuebiq provided mobility data for a number of COVID-19 studies that summarize the distance users travel or the proportion of users that stay at home52. These existing frameworks – not only the analyses, but also the privacy considerations and data sharing agreements – will provide standardized methods that facilitate integrating mobility data into intervention assessments.
Data privacy
Various forms of identifiable personal information are generated when using mobile phones, including names, identification numbers, fine spatial and temporal data on where the device was used, other users’ identification numbers who may have been detected by Bluetooth, and personal details that might be entered into an app. In light of the growing number of digital privacy concerns and regulations, one must carefully consider the exact form and use of mobile phone data being collected against the legal and ethical need to protect users’ data security and confidentiality. While maintaining user confidentiality is often seen as a hindrance to the use of mobile phone data, in that it limits the use of individual-level data and typically requires aggregation to coarse spatial and temporal resolutions, there are a number of existing frameworks that can help provide guidance for the effective, privacy-conscious use of mobile phone data53.
Exactly which model of data privacy will best suit the use of mobile phone data for COVID-19 response will depend on the exact form and proposed use of the data. As discussed above, there already exist many data processing and analysis frameworks to provide anonymized indicators of population mobility. These standard procedures, though, could result in aggregated data with insufficient spatial and temporal resolution to be effective for monitoring the spread of SARS-CoV-2. Privacy regulations, such as the European Union’s General Data Protection Regulation (GDPR)54, offer exceptions for the use of non-anonymous data that may be needed for other response efforts. For example, opt-in applications for contact tracing may seek consent of the data subject to collect and analyze identifiable data, though the ability to scale opt-in approaches to a wide enough population and to maintain user compliance and participation remains unclear. GDPR and other regulations also provide an exception for anonymization of data to be used in public service, but the regulatory hurdles to gain this exception can be substantial and would require clear use policies and applications for these data. The use of mobile phone data, particularly forms such as those proposed through contact tracing applications, must be weighed against the possible infringements of privacy and civil liberties versus the potential public health benefit. More

Land use mapping and accuracy assessment
According to the land use planning map of Zhuhai city, the characteristics of the city, the status of human activities and land use, and the types of natural ecosystems, we identified and categorized land use into 10 types: woodland, grassland, rainfed cropland, paddy fields, aquaculture areas, reservoirs and pit ponds, tidal flats, rivers and shallow water, built-up land and unutilized land (Supplemental Materials S1: Land use types and descriptions). The ecological land types include woodland, grassland, reservoirs and pit ponds, tidal flats, and rivers and shallow water. Rainfed cropland, paddy fields and aquaculture areas were not included as ecological land types because they are agricultural land mainly used for agricultural production. These land use types are greatly disturbed by humans, their ecological functions are very fragile, and they are affected by economic interests and have low ecological value. Unutilized land provides few ecological benefits and may be converted into built-up land in the short term; thus, its ecological benefits are unsustainable.
After the preprocessing and splicing of multiperiod satellite RS images, we completed object-based multiscale automatic segmentation and land use classification of the images using eCognition Developer software. Specifically, the Estimation of Scale Parameters (ESP) tool was first used to obtain the local variance parameter, which reflects the internal homogeneity of the segmentation object; then, the rate of change (ROC) of the local variance (LV) parameter was calculated37,38. When the ROC reaches its peak, the corresponding segmentation scale can be used as the optimal segmentation scale37. At the optimal segmentation scale, classification is based on the object unit using the nearest neighbor method of eCognition Developer. The nearest neighbor method is a commonly used supervised classification method that is simple and easy to understand, and it is suitable for multiclassification problems39.
Finally, based on the preliminary results data of the four stages automatically classified by eCognition Developer, obvious errors and omissions in the data of the preliminary results were revised and improved through manual visual interpretation. The final revised result data were used for the subsequent analysis of the land pattern and its changes.
This study first drew land use maps for four years: 1991, 2000, 2010, and 2018. We extracted no less than 200 regions of interest (ROIs) in each study year and compared high-resolution Google Earth images to perform a land mapping accuracy assessment. To ensure that the accuracy of each land type was reliably estimated, we confirmed that each land type had at least 10 ROIs when laying out the ROI area. Table 1 shows the land use classification accuracy for the 1991–2018 period. The overall accuracy of the land mapping for 1991, 2000, 2010, and 2018 was 93.4%, 94.1%, 91.1%, and 94.5%, respectively, and the Kappa coefficients were 0.925, 0.933, 0.890, and 0.938, respectively, meeting the research requirements.
Table 1 Classification accuracy of land use types in Zhuhai city.
Full size table

Spatial patterns and dynamics of ecological land
From 1991 to 2018, the ecological land in Zhuhai was dominated by woodland and rivers and shallow water, and the overall area of ecological land continuously decreased (Fig. 1). In 1991, the total area of ecological land was 849.4 km2, accounting for 53.7% of Zhuhai’s urban area. In 2018, the area was reduced to 574.6 km2, accounting for only 36.3% of Zhuhai’s urban area.
Figure 1

The net change in ecological land in Zhuhai city, 1991–2018. The area of woodland is the largest, followed by the area of rivers and shallow water. The proportions of woodland and grassland in the total area of ecological land increased by 7.6% and 1.3%, respectively. Rivers and shallow water and tidal flats showed downward trends, decreasing by 8.7% and 1.8%, respectively. Reservoirs and pit ponds increased slightly and showed dynamic changes.

Full size image

In 28 years, the amount of ecological land decreased by 32.3%, of which woodland decreased by 24.2% (129.6 km2), tidal flats decreased by 67.2% (19.3 km2), and rivers and shallow water decreased by as much as 51.8% (132.3 km2). The reduction in rivers and shallow water represented the bulk of the reduction in ecological land area (48.1%). In contrast, the area of reservoirs and pit ponds grew slightly while maintaining a steady state, increasing by 1.1 km2. Compared with 1991, the grassland area grew slightly, increasing by 5.3 km2, mainly due to the construction of golf courses and parks. Clearly, there is an order of magnitude difference between the increase and decrease in ecological land.
From the temporal perspective (Fig. 2), the change in ecological land mainly occurred in the 1991–2000 period. During this period, the reduction in ecological land was the largest (212.3 km2), mainly distributed in the contiguous area of woodland and built-up land in the central and western areas of the Doumen District and in the coastal areas of the Jinwan District and Xiangzhou District. At the same time, there was a small increase in ecological land, mainly due to the restoration and regulation of tidal flats and reservoirs and pit ponds.
Figure 2

Ecological land gains and losses in Zhuhai city, 1991–2018. (a,c,e) show an increase in ecological land; (b,d,f) show a decrease in ecological land. The decrease in ecological land is obviously higher than the increase, and there is an increase in the degree of patch fragmentation. The reduced patches are mostly marginal woodland and river and shallow water areas. The boundaries of the map come from the Zhuhai Natural Resources Bureau. The drawing of the map was completed with the support of ArcGIS 10.7 software.

Full size image

Since 2000, ecological environmental protection and construction work have gradually been taken more seriously, and the State Council of China promulgated the “National Ecological Environmental Protection Program”. Local governments at all levels have gradually strengthened their awareness of ecological environmental protection. The occupation of ecological land by urban development has rapidly decreased, while the area of new ecological land formed by ecological protection and ecological restoration has gradually and steadily increased. From 2000 to 2010, the ecological land in Zhuhai decreased by 130.1 km2 and increased by 53.6 km2, with a net reduction of 76.5 km2. From 2010 to 2018, the decrease and increase in ecological land were similar, and the net reduction in area was only 18.6 km2; thus, the spatial distribution and quantity of ecological land in Zhuhai city was approximately stable (Fig. 3).
Figure 3

Losses and gains in ecological land area in Zhuhai city, 1991–2018. Green indicates an increase in ecological land, and red indicates a decrease in ecological land. From 1991 to 2000, the net reduction in ecological land was 177.9 km2. From 2000 to 2010, the net reduction in ecological land was 76.5 km2. From 2010 to 2018, the net reduction in ecological land was 18.6 km2.

Full size image

In the 28-year monitoring period of this paper, the reduction in ecological land in the first 10 years (1991–2000) was 0.99 times that in the subsequent 18 years (2000–2018). The total amount of ecological land added in the subsequent 18 years (2000–2018) was 3.5 times that of the first 10 years (1991–2000).
Landscape characteristics
At the landscape level (Table 2), the edge density (ED) of ecological land in the study area is significantly lower than that of nonecological land. The ED exhibited a pattern of first increasing, then decreasing, and subsequently slightly increasing (with values of 33.6 in 1991, 37.7 in 2000, 31.8 in 2010, and 34.7 in 2018). The patch density (PD), landscape shape index (LSI), and largest patch index (LPI) had the same trend as that of the ED. These changes indicate that over time, the landscape of ecological land began to experience an increase in fragmentation and a decrease in regularity and continuity; then, the landscape was reintegrated into a more regular and continuous pattern.
Table 2 Changes in landscape-level indexes in Zhuhai city, 1991–2018.
Full size table

In addition, from 1991 to 2018, the contagion index (CONTAG) of all land in Zhuhai city fluctuated slightly at approximately 55%, and the degree of landscape pattern aggregation did not change much. However, the CONTAG of ecological land was approximately 70%, which was significantly higher than that of nonecological land; this result indicates that the CONTAG and connectivity of ecological land were higher than those of nonecological land. Shannon’s diversity index (SHDI) and Shannon’s evenness index (SHEI) did not change much in the time series, indicating that the landscape diversity of Zhuhai city has basically been stable over the past 28 years. However, compared with 1991, the SHDI and SHEI decreased slightly, indicating that the ecological landscape diversity and uniformity decreased in the study area, while the landscape heterogeneity increased.
At the class level (Table 3), the PD and the area-weighted mean contiguity index (CONTIG_AM) of woodland remained basically unchanged, the LSI increased from 19.99 to 21.7, and the LPI decreased from 9.6 to 3.9. These changes were caused by the following processes: the expansion of built-up land, the preferential occupation of marginal forestland by built-up land, the reduction in the dominance of the landscape type, and the increasing complexity of the original geometry. However, woodland mainly exists in a continuous form, and these encroachment behaviors have little effect on the number, spatial connectivity or proximity of woodland patches.
Table 3 Changes in class-level indexes in Zhuhai city, 1991–2018.
Full size table

The PD and LSI of grassland showed downward trends, while the LPI and CONTIG_AM showed upward trends. This result is closely related to the increase in grassland in the study area. The increased grassland caused the number of patches to increase slightly, improving the superiority of the landscape. The construction of artificial grassland is more regular in the shape of grass patches, and the connectivity is enhanced between landscape units.
In addition, the PD, LSI and LPI of tidal flats showed downward trends, indicating that the development and utilization of tidal flat reclamation were strengthened, the number decreased, and the shape tended to be regular. The landscape characteristics of reservoirs and pit ponds and rivers and shallow water were basically the same: the LSI showed an upward trend, indicating that the patches were seriously disturbed by human activities, the large patches experienced continuous fragmentation, and the landscape type shapes were complicated. In contrast, the LPI showed a downward trend, indicating that activities such as sea filling led to a continuous decrease in sea area.
Ecological quality evaluation
Ecological quality is used to characterize the conditions of the ecosystem; the ecosystem is disturbed by human activities and land use change, and the ability to provide services is also affected40. The value of ecosystem services is an important comprehensive indicator reflecting ecological quality, and the ecological service value of ecological land is higher than that of nonecological land41. Based on the ecosystem service value coefficient proposed by Xie et al.28, we normalized the coefficient value to 0–1 and used the equivalent area and the average equivalent area, which were used to evaluate the ecological service quality of ecological land.
The transformation matrix of ecological land and nonecological land shows the following (Table 4): the probability of ecological land being transformed into nonecological land in the periods 1991–2000, 2000–2010 and 2010–2018 was 25.0%, 19.4% and 14.3%, respectively. The contributions of ecological land to nonecological land were 23.3%, 13.2% and 8.4%, respectively. The transformation of ecological land to nonecological land showed a weakening trend after 2000, and the ecological quality showed improvement.
Table 4 Probability of ecological land being transformed into nonecological land in Zhuhai city, 1991–2018.
Full size table

From 1991 to 2018, the equivalent area of ecological land continued to decrease, but the downward trend gradually stabilized after 2000 (Fig. 4). In 1991, the equivalent area of regional ecological land was 849.4 km2, and in 2000, it was 673.2 km2, indicating a significant decrease in the equivalent area, with a reduction of 20.7%. In 2010, the equivalent area of ecological land further dropped to 600.2 km2, a reduction of 10.8%, although the decrease was significantly smaller than that in the previous period. In 2018, the equivalent area was 574.6 km2, representing a reduction of only 4.3%.
Figure 4

Dynamic changes in ecological land quality in Zhuhai city, 1991–2018. From 1991 to 2018, the equivalent area of ecological land in Zhuhai city showed a downward trend, with a decrease of 274.8 km2, i.e., 32.3%. The average equivalent area index showed an upward trend, with an increase of 0.11, i.e., 9.3%.

Full size image

As shown in Fig. 4, the average equivalent area of ecological land showed a continuous upward trend. Specifically, the average equivalent area was 1.14 in 1991, 1.22 in 2000, 1.24 in 2010, and 1.25 in 2018. This result shows that although the ecological land area decreased, the quality of the ecological land gradually improved. In reality, this pattern was manifested as follows: the area of grasslands and reservoirs and pit ponds gradually increased, the degree of landscape fragmentation weakened, and the landscape dominance became more obvious. In addition, these land types have relatively high ecosystem service values among all land types.
Changes in the center of gravity of ecological land
From 1991 to 2018, the center of gravity of ecological land shifted to the northeast, and the center of gravity of built-up land shifted to the southwest (Fig. 5).
Figure 5

Changes in the center of gravity of ecological land and built-up land in Zhuhai city, 1991–2018. From 1991 to 2018, the center of gravity of ecological land in Zhuhai moved to the northeast by 1346 m. The center of gravity of built-up land moved in the opposite direction, moving 7254 m to the southwest. The boundaries of the map come from the Zhuhai Natural Resources Bureau, and the base map in the main map is the China Online Community Basemap in ArcGIS. The drawing of this map was completed with the support of ArcGIS 10.7 software.

Full size image

From 1991 to 2000, the center of gravity of ecological land moved 404 m to the east and 409 m to the north, and the overall movement was 578 m to the northeast. From 2000 to 2010, the center of gravity of ecological land moved 24 m to the east and 355 m to the north, and the overall movement trend was northward. From 2010 to 2018, the center of gravity of ecological land moved 273 m to the east and 236 m to the north, and the overall movement was 473 m to the northeast. In these three periods, the center of gravity of built-up land moved to the southwest by 2871 m, 3983 m and 424 m. The urban expansion and internal construction mainly experienced a rapid and then slow evolution from the northeast to the southwest.
From the spatial distribution of all ecological land types, the center of gravity of woodland moved to the southeast (0.68 km) from 1991 to 2018. This movement occurred because human construction activities such as deforestation, urban expansion, and infrastructure construction were prominent in the western and northern parts of Zhuhai during the 1991–2000 period. The movement of the center of gravity of grassland to the east and south was highly related to the construction of golf courses, such as the Zhuxiandong Golf Club in the Xiangzhou District, the Dananshan Cuihu Golf Course in Jinding Town, a golf club in the Jinwan District, and Zhuhai Stadium in the Xiangzhou District. The center of gravity of reservoirs and pit ponds moved southward (2.9 km); the center of gravity of tidal flats moved eastward (5.8 km); and the center of gravity of rivers and shallow water moved northward (3.5 km). These changes were closely related to the reclamation engineering carried out by Zhuhai city in recent years.
Modeling the ecological land change process
Changes in urban ecological land are mainly due to the expansion of the outer edge of cities and the oppression of urban internal land development. Therefore, we selected four indicators of natural geography and regional development that might reflect changes in urban expansion and urban construction: elevation, slope, distance from built-up land, and growth rate of built-up land.
With the support of SPSS software, the equation of the transformation probability of ecological land to nonecological land in Zhuhai can be obtained through the binary logistic regression analysis module. Specifically, this equation is expressed as follows (see Supplemental Materials S2: Parameter of the driving factors for modeling):

$$P = 1 – frac{1}{{{1 + }e^{{{ – }left( {{0}{text{.069}} times {text{A } + text{ 0}}{.033} times {text{B } + text{ 0}}{.473} times {text{C } – text{ 1}}{.079} times {text{D } – text{ 0}}{.963} times {text{E } – text{ 0}}{.853}} right)}} }}$$
(1)

where A is the slope; B is the elevation; C is the distance from built-up land; and D and E are the built-up land growth rates of categories 4 and 5, respectively. The squared maximum likelihood of the numerical values (− 2 log-likelihood) of the model was 18,155.4, and the value of the χ2(5) comprehensive test statistic was 7871.2 (p  More

1.
Nowell RW, Green S, Laue BE, Sharp PM. The extent of genome flux and its role in the differentiation of bacterial lineages. Genome Biol Evol. 2014;6:1514–29.
PubMed  PubMed Central  Article  CAS  Google Scholar 
2.
Ochman H, Lawrence JG, Groisman EA. Lateral gene transfer and the nature of bacterial innovation. Nature. 2000;405:299–304.
CAS  PubMed  Article  Google Scholar 

3.
Wiedenbeck J, Cohan FM. Origins of bacterial diversity through horizontal genetic transfer and adaptation to new ecological niches. FEMS Microbiol Rev. 2011;35:957–76.
CAS  PubMed  PubMed Central  Article  Google Scholar 

4.
Albalat R, Cañestro C. Evolution by gene loss. Nat Rev Genet. 2016;17:379–91.
CAS  PubMed  Article  Google Scholar 

5.
Jacq C, Miller JR, Brownlee GG. A pseudogene structure in 5S DNA of Xenopus laevis. Cell. 1977;12:109–20.
CAS  PubMed  Article  Google Scholar 

6.
Li W-H, Gojobori T, Nei M. Pseudogenes as a paradigm of neutral evolution. Nature. 1981;292:237–9.
CAS  PubMed  Article  Google Scholar 

7.
Ohta T. The nearly neutral theory of molecular evolution. Annu Rev Ecol Syst. 1992;23:263–86.
Article  Google Scholar 

8.
Bolotin E, Hershberg R. Gene loss dominates as a source of genetic variation within clonal pathogenic bacterial species. Genome Biol Evol. 2015;7:2173–87.
CAS  PubMed  PubMed Central  Article  Google Scholar 

9.
Hottes AK, Freddolino PL, Khare A, Donnell ZN, Liu JC, Tavazoie S. Bacterial adaptation through loss of function. PLoS Genet. 2013;9:e1003617.
CAS  PubMed  PubMed Central  Article  Google Scholar 

10.
Sharma V, Hecker N, Roscito JG, Foerster L, Langer BE, Hiller M. A genomics approach reveals insights into the importance of gene losses for mammalian adaptations. Nat Commun. 2018;9:1–9.
Article  CAS  Google Scholar 

11.
Sokurenko EV, Hasty DL, Dykhuizen DE. Pathoadaptive mutations: gene loss and variation in bacterial pathogens. Trends Microbiol. 1999;7:191–5.
CAS  PubMed  Article  Google Scholar 

12.
Ortega AP, Villagra NA, Urrutia IM, Valenzuela LM, Talamilla-Espinoza A, Hidalgo AA, et al. Lose to win: marT pseudogenization in Salmonella enterica serovar Typhi contributed to the surV-dependent survival to H2O2, and inside human macrophage-like cells. Infect Genet Evol. 2016;45:111–21.
CAS  PubMed  Article  Google Scholar 

13.
Goodhead I, Darby AC. Taking the pseudo out of pseudogenes. Curr Opin Microbiol. 2015;23:102–9.
CAS  PubMed  Article  Google Scholar 

14.
Johnson LJ. Pseudogene rescue: an adaptive mechanism of codon reassignment. J Evol Biol. 2010;23:1623–30.
CAS  PubMed  Article  Google Scholar 

15.
Librado P, Vieira FG, Rozas J. BadiRate: estimating family turnover rates by likelihood-based methods. Bioinformatics. 2012;28:279–81.
CAS  PubMed  Article  Google Scholar 

16.
David LA, Alm EJ. Rapid evolutionary innovation during an Archaean genetic expansion. Nature. 2011;469:93–96.
CAS  PubMed  Article  Google Scholar 

17.
Avni E, Montoya D, Lopez D, Modlin R, Pellegrini M, Snir S. A phylogenomic study quantifies competing mechanisms for pseudogenization in prokaryotes—the Mycobacterium leprae case. PLoS One. 2017;13:e0204322.
Article  CAS  Google Scholar 

18.
Ochman H. The nature and dynamics of bacterial genomes. Science. 2006;311:1730–3.
CAS  PubMed  Article  Google Scholar 

19.
Grote J, Thrash JC, Huggett MJ, Landry ZC, Carini P, Giovannoni SJ, et al. Streamlining and core genome conservation among highly divergent members of the SAR11 clade. mBio. 2012;3:e00252–12.
CAS  PubMed  PubMed Central  Article  Google Scholar 

20.
Giovannoni SJ, Cameron Thrash J, Temperton B. Implications of streamlining theory for microbial ecology. ISME J. 2014;8:1553–65.
PubMed  PubMed Central  Article  Google Scholar 

21.
Buchan A, González JM, Moran MA. Overview of the marine Roseobacter lineage. Appl Environ Microbiol. 2005;71:5665–77.
CAS  PubMed  PubMed Central  Article  Google Scholar 

22.
Luo H, Moran MA. How do divergent ecological strategies emerge among marine bacterioplankton lineages? Trends Microbiol. 2015;23:577–84.
CAS  PubMed  Article  Google Scholar 

23.
Luo H, Moran MA. Evolutionary ecology of the marine Roseobacter clade. Microbiol Mol Biol Rev. 2014;78:573–87.
PubMed  PubMed Central  Article  Google Scholar 

24.
Tujula NA, Crocetti GR, Burke C, Thomas T, Holmström C, Kjelleberg S. Variability and abundance of the epiphytic bacterial community associated with a green marine Ulvacean alga. ISME J. 2010;4:301–11.
PubMed  Article  Google Scholar 

25.
Littman RA, Willis BL, Pfeffer C, Bourne DG. Diversities of coral-associated bacteria differ with location, but not species, for three acroporid corals on the Great Barrier Reef. FEMS Microbiol Ecol. 2009;68:152–63.
CAS  PubMed  Article  Google Scholar 

26.
Rosenberg E, Koren O, Reshef L, Efrony R, Zilber-Rosenberg I. The role of microorganisms in coral health, disease and evolution. Nat Rev Microbiol. 2007;5:355–62.
CAS  PubMed  Article  Google Scholar 

27.
Sweet MJ, Croquer A, Bythell JC. Bacterial assemblages differ between compartments within the coral holobiont. Coral Reefs. 2011;30:39–52.
Article  Google Scholar 

28.
Crossland CJ, Barnes DJ, Borowitzka MA. Diurnal lipid and mucus production in the staghorn coral Acropora acuminata. Mar Biol. 1980;60:81–90.
CAS  Article  Google Scholar 

29.
Shashar N, Stambler N. Endolithic algae within corals—life in an extreme environment. J Exp Mar Biol Ecol. 1992;163:277–86.
CAS  Article  Google Scholar 

30.
Highsmith RC. Lime-boring algae in hermatypic coral skeletons. J Exp Mar Biol Ecol. 1981;55:267–81.
Article  Google Scholar 

31.
Kühl M, Holst G, Larkum AWD, Ralph PJ. Imaging of oxygen dynamics within the endolithic algal community of the massive coral Porites lobata. J Phycol. 2008;44:541–50.
PubMed  Article  CAS  Google Scholar 

32.
Kalhoefer D, Thole S, Voget S, Lehmann R, Liesegang H, Wollher A, et al. Comparative genome analysis and genome-guided physiological analysis of Roseobacter litoralis. BMC Genomics. 2011;12:324.
CAS  PubMed  PubMed Central  Article  Google Scholar 

33.
Lachnit T, Fischer M, Künzel S, Baines JF, Harder T. Compounds associated with algal surfaces mediate epiphytic colonization of the marine macroalga Fucus vesiculosus. FEMS Microbiol Ecol. 2013;84:411–20.
CAS  PubMed  Article  Google Scholar 

34.
Singh RP, Reddy CRK. Seaweed–microbial interactions: key functions of seaweed-associated bacteria. FEMS Microbiol Ecol. 2014;88:213–30.
CAS  PubMed  Article  Google Scholar 

35.
Khailov KM, Burlakova ZP. Release of dissolved organic matter by marine seaweeds and distribution of their total organic production to inshore communities. Limnol Oceanogr. 1969;14:521–7.
Article  Google Scholar 

36.
Wai TC, Ng JSS, Leung KMY, Dudgeon D, Williams GA. The source and fate of organic matter and the significance of detrital pathways in a tropical coastal ecosystem. Limnol Oceanogr. 2008;53:1479–92.
CAS  Article  Google Scholar 

37.
Braeckman U, Pasotti F, Vázquez S, Zacher K, Hoffmann R, Elvert M, et al. Degradation of macroalgal detritus in shallow coastal Antarctic sediments. Limnol Oceanogr. 2019;64:1423–41.
CAS  PubMed  PubMed Central  Google Scholar 

38.
Moran MA, Belas R, Schell MA, Gonzalez JM, Sun F, Sun S, et al. Ecological genomics of marine roseobacters. Appl Environ Microbiol. 2007;73:4559–69.
CAS  PubMed  PubMed Central  Article  Google Scholar 

39.
Sonnenschein EC, Nielsen KF, D’Alvise P, Porsby CH, Melchiorsen J, Heilmann J, et al. Global occurrence and heterogeneity of the Roseobacter clade species Ruegeria mobilis. ISME J. 2017;11:569–83.
CAS  PubMed  Article  PubMed Central  Google Scholar 

40.
Slightom RN, Buchan A. Surface colonization by marine roseobacters: integrating genotype and phenotype. Appl Environ Microbiol. 2009;75:6027–37.
CAS  PubMed  PubMed Central  Article  Google Scholar 

41.
Thole S, Kalhoefer D, Voget S, Berger M, Engelhardt T, Liesegang H, et al. Phaeobacter gallaeciensis genomes from globally opposite locations reveal high similarity of adaptation to surface life. ISME J. 2012;6:2229–44.
CAS  PubMed  PubMed Central  Article  Google Scholar 

42.
Newton RJ, Griffin LE, Bowles KM, Meile C, Gifford S, Givens CE, et al. Genome characteristics of a generalist marine bacterial lineage. ISME J. 2010;4:784–98.
CAS  PubMed  Article  Google Scholar 

43.
Brinkhoff T, Giebel H-A, Simon M. Diversity, ecology, and genomics of the Roseobacter clade: a short overview. Arch Microbiol. 2008;189:531–9.
CAS  PubMed  Article  Google Scholar 

44.
Luo H, Löytynoja A, Moran MA. Genome content of uncultivated marine Roseobacters in the surface ocean. Environ Microbiol. 2012;14:41–51.
CAS  PubMed  Article  Google Scholar 

45.
Lerat E, Ochman H. Ψ-Φ: Exploring the outer limits of bacterial pseudogenes. Genome Res. 2004;14:2273–8.
CAS  PubMed  PubMed Central  Article  Google Scholar 

46.
Lerat E, Ochman H. Recognizing the pseudogenes in bacterial genomes. Nucleic Acids Res. 2005;33:3125–32.
CAS  PubMed  PubMed Central  Article  Google Scholar 

47.
Kuo C-H, Ochman H. The extinction dynamics of bacterial pseudogenes. PLoS Genet. 2010;6:e1001050.
PubMed  PubMed Central  Article  CAS  Google Scholar 

48.
Umezaki I. Ecological studies of Sargassum hemiphyllum C. AGARDH in Obama Bay, Japan Sea. Nippon Suisan Gakkaishi. 1984;50:1677–83.
Article  Google Scholar 

49.
Tam TW, Ang PO. Repeated physical disturbances and the stability of sub-tropical coral communities in Hong Kong. China Aquat Conserv Mar Freshw Ecosyst. 2008;18:1005–24.
Article  Google Scholar 

50.
Cheang CC, Chu KH, Ang PO. Phylogeography of the marine macroalga Sargassum hemiphyllum (Phaeophyceae, Heterokontophyta) in northwestern Pacific. Mol Ecol. 2010;19:2933–48.
CAS  PubMed  Article  Google Scholar 

51.
Raghunathan C, Venkataraman K. Diversity and distribution of corals and their associated fauna of Rani Jhansi Marine National Park, Andaman and Nicobar Islands. In: Venkataraman K, Raghunathan C, Sivaperuman C, (eds). Ecology of Faunal Communities on the Andaman and Nicobar Islands. Berlin, Heidelberg: Springer; 2012. p. 177–208.
Google Scholar 

52.
Ang PO. Phenology of Sargassum spp. in Tung Ping Chau Marine Park, Hong Kong SAR, China. J Appl Phycol. 2006;18:629–36.
Article  Google Scholar 

53.
Huggett MJ, Apprill A. Coral microbiome database: Integration of sequences reveals high diversity and relatedness of coral-associated microbes. Environ Microbiol Rep. 2019;11:372–85.
PubMed  Article  Google Scholar 

54.
Passel MWJ, van, Marri PR, Ochman H. The emergence and fate of horizontally acquired genes in Escherichia coli. PLoS Comput Biol. 2008;4:e1000059.
PubMed  PubMed Central  Article  CAS  Google Scholar 

55.
Ochman H. Distinguishing the ORFs from the ELFs: short bacterial genes and the annotation of genomes. Trends Genet. 2002;18:335–7.
CAS  PubMed  Article  Google Scholar 

56.
Emms DM, Kelly S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 2019;20:238.
PubMed  PubMed Central  Article  Google Scholar 

57.
Nguyen LT, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32:268–74.
CAS  Article  Google Scholar 

58.
Schliep KP. Phangorn: phylogenetic analysis in R. Bioinformatics. 2011;27:592–3.
CAS  PubMed  PubMed Central  Article  Google Scholar 

59.
Liu Y, Harrison PM, Kunin V, Gerstein M. Comprehensive analysis of pseudogenes in prokaryotes: widespread gene decay and failure of putative horizontally transferred genes. Genome Biol. 2004;5:R64.
PubMed  PubMed Central  Article  Google Scholar 

60.
Halldal P. Photosynthetic capacities and photosynthetic action spectra of endozoic algae of the massive coral Favia. Biol Bull. 1968;134:411–24.
CAS  Article  Google Scholar 

61.
Shibata K, Haxo FT. Light transmission and spectral distribution through epi- and endozoic algal layers in the brain coral, Favia. Biol Bull. 1969;136:461–8.
CAS  Article  Google Scholar 

62.
Park JT, Uehara T. How bacteria consume their own exoskeletons (turnover and recycling of cell wall peptidoglycan). Microbiol Mol Biol Rev. 2008;72:211–27.
CAS  PubMed  PubMed Central  Article  Google Scholar 

63.
Mauck J, Chan L, Glaser L. Turnover of the cell wall of gram-positive bacteria. J Biol Chem. 1971;246:1820–7.
CAS  PubMed  Google Scholar 

64.
Goodell E. Recycling of murein by Escherichia coli. J Bacteriol. 1985;163:305–10.
CAS  PubMed  PubMed Central  Article  Google Scholar 

65.
Uehara T, Suefuji K, Jaeger T, Mayer C, Park JT. MurQ etherase is required by Escherichia coli in order to metabolize Anhydro-N-Acetylmuramic acid obtained either from the environment or from its own cell wall. J Bacteriol. 2006;188:1660–2.
CAS  PubMed  PubMed Central  Article  Google Scholar 

66.
Dik DA, Marous DR, Fisher JF, Mobashery S. Lytic transglycosylases: concinnity in concision of the bacterial cell wall. Crit Rev Biochem Mol Biol. 2017;52:503–42.
PubMed  PubMed Central  Article  Google Scholar 

67.
Jiang H, Kong R, Xu X. The N-acetylmuramic acid 6-phosphate etherase gene promotes growth and cell differentiation of cyanobacteria under light-limiting conditions. J Bacteriol. 2010;192:2239–45.
CAS  PubMed  PubMed Central  Article  Google Scholar 

68.
Ferrer LM, Szmant AM. Nutrient regeneration by the endolithic community in coral skeletons. In: Proceedings of the 6th International Coral Reef Symposium 1988. pp 1–4.

69.
Risk MJ, Muller HR. Porewater in coral heads: evidence for nutrient regeneration. Limnol Oceanogr. 1983;28:1004–8.
Article  Google Scholar 

70.
Yu LJ, Wu JR, Zheng ZZ, Lin CC, Zhan XB. Changes in gene transcription and protein expression involved in the response of Agrobacterium sp. ATCC 31749 to nitrogen availability during curdlan production. Appl Biochem Microbiol. 2011;47:487–93.
CAS  Article  Google Scholar 

71.
Wada S, Aoki M, Mikami A, Komatsu T, Tsuchiya Y, Sato T, et al. Bioavailability of macroalgal dissolved organic matter in seawater. Mar Ecol Prog Ser. 2008;370:33–44.
CAS  Article  Google Scholar 

72.
Essenberg MK, Cooper RA. Two ribose-5-phosphate isomerases from Escherichia coli K12: partial characterisation of the enzymes and consideration of their possible physiological roles. Eur J Biochem. 1975;55:323–32.
CAS  PubMed  Article  Google Scholar 

73.
Nelson CE, Goldberg SJ, Wegley Kelly L, Haas AF, Smith JE, Rohwer F, et al. Coral and macroalgal exudates vary in neutral sugar composition and differentially enrich reef bacterioplankton lineages. ISME J. 2013;7:962–79.
CAS  PubMed  PubMed Central  Article  Google Scholar 

74.
Mulligan C, Fischer M, Thomas GH. Tripartite ATP-independent periplasmic (TRAP) transporters in bacteria and archaea. FEMS Microbiol Rev. 2011;35:68–86.
CAS  PubMed  Article  Google Scholar 

75.
Beyenbach KW, Wieczorek H. The V-type H+ ATPase: molecular structure and function, physiological roles and regulation. J Exp Biol. 2006;209:577–89.
CAS  PubMed  Article  Google Scholar 

76.
Guadayol Ò, Silbiger NJ, Donahue MJ, Thomas FIM. Patterns in temporal variability of temperature, oxygen and pH along an environmental gradient in a coral reef. PLoS One. 2014;9:e85213.
PubMed  PubMed Central  Article  CAS  Google Scholar 

77.
Bodenmiller DM, Spiro S. The yjeB(nsrR) gene of Escherichia coli encodes a nitric oxide-sensitive transcriptional regulator. J Bacteriol. 2006;188:874–81.
CAS  PubMed  PubMed Central  Article  Google Scholar 

78.
Gilberthorpe NJ, Lee ME, Stevanin TM, Read RC, Poole RK. NsrR: a key regulator circumventing Salmonella enterica serovar Typhimurium oxidative and nitrosative stress in vitro and in IFN-γ-stimulated J774.2 macrophages. Microbiology. 2007;153:1756–71.
CAS  PubMed  PubMed Central  Article  Google Scholar 

79.
da Fonseca RR, Johnson WE, O’Brien SJ, Vasconcelos V, Antunes A. Molecular evolution and the role of oxidative stress in the expansion and functional diversification of cytosolic glutathione transferases. BMC Evol Biol. 2010;10:281.
PubMed  PubMed Central  Article  CAS  Google Scholar 

80.
Green ER, Mecsas J. Bacterial secretion systems—an overview. Microbiol Spectr. 2016;4:1–32.
CAS  Article  Google Scholar 

81.
Ansari MI, Schiwon K, Malik A, Grohmann E. Biofilm formation by environmental bacteria. In: Malik A, Grohmann E (eds). Environmental protection strategies for sustainable development. 2012. Springer Netherlands, Dordrecht, pp 341–77.

82.
Meron D, Efrony R, Johnson WR, Schaefer AL, Morris PJ, Rosenberg E, et al. Role of flagella in virulence of the coral pathogen Vibrio coralliilyticus. Appl Environ Microbiol. 2009;75:5704–7.
CAS  PubMed  PubMed Central  Article  Google Scholar 

83.
Attmannspacher U, Scharf BE, Harshey RM. FliL is essential for swarming: motor rotation in absence of FliL fractures the flagellar rod in swarmer cells of Salmonella enterica. Mol Microbiol. 2008;68:328–41.
CAS  PubMed  Article  Google Scholar 

84.
Fernando SC, Wang J, Sparling K, Garcia GD, Francini-Filho RB, de Moura RL, et al. Microbiota of the major south atlantic reef building coral Mussismilia. Micro Ecol. 2015;69:267–80.
Article  Google Scholar 

85.
Pollock FJ, McMinds R, Smith S, Bourne DG, Willis BL, Medina M, et al. Coral-associated bacteria demonstrate phylosymbiosis and cophylogeny. Nat Commun. 2018;9:4921.
PubMed  PubMed Central  Article  CAS  Google Scholar 

86.
Marcelino VR, van Oppen MJ, Verbruggen H. Highly structured prokaryote communities exist within the skeleton of coral colonies. ISME J. 2018;12:300–3.
PubMed  Article  Google Scholar 

87.
Hill C. Virulence or niche factors: what’s in a name? J Bacteriol. 2012;194:5725–7.
CAS  PubMed  PubMed Central  Article  Google Scholar 

88.
Egan S, Harder T, Burke C, Steinberg P, Kjelleberg S, Thomas T. The seaweed holobiont: understanding seaweed–bacteria interactions. FEMS Microbiol Rev. 2013;37:462–76.
CAS  PubMed  Article  Google Scholar 

89.
Levy A, Salas Gonzalez I, Mittelviefhaus M, Clingenpeel S, Herrera Paredes S, Miao J, et al. Genomic features of bacterial adaptation to plants. Nat Genet. 2018;50:138–50.
CAS  Article  Google Scholar 

90.
Koren O, Rosenberg E. Bacteria associated with mucus and tissues of the coral Oculina patagonica in summer and winter. Appl Environ Microbiol. 2006;72:5254–9.
CAS  PubMed  PubMed Central  Article  Google Scholar 

91.
Wang X, Grus WE, Zhang J. Gene losses during human origins. PLoS Biol. 2006;4:e52.
PubMed  PubMed Central  Article  CAS  Google Scholar  More

Inhibition of cell signaling by commensal streptococci
To study how S. mutans ComRS signaling could be impacted by the presence of a competing species, we empirically optimized a dual-species model system (Fig. 1a) in which a strain of S. mutans carrying the promoter regions of comS or comX (PcomS, PcomX) fused to a codon-optimized green fluorescence protein (gfp) reporter gene could be cocultured with wild-type strains of Streptococcus gordonii DL1, Streptococcus sanguinis SK150, or S. sp. A12. All experiments were performed in chemically defined medium (CDM) [38, 39] because activation of the ComRS circuit occurs spontaneously in CDM as cell density increases, with no need for addition of synthetic XIP or overexpression of the gene for the XIP precursor (comS) (Supplementary Fig. 1). CDM is also heavily buffered with phosphate, which is advantageous because ComRS signaling is optimal at neutral pH values [40, 41]. The buffer also prevents the generation of strongly acidic conditions by S. mutans, which is detrimental to the comparatively acid-sensitive commensal Streptococcus spp.
Fig. 1: Loss of S. mutans peptide signaling in presence of competitor.

a An oral Streptococcus spp. competitor strain (blue) was cocultured in chemically defined medium (CDM) with an S. mutans PcomX::gfp reporter strain (green). As cell density of the reporter strain increases during growth, the XIP peptide that originates from the comS gene will be produced and accumulates extracellularly. XIP is then reimported into the cell through the Opp oligopeptide permease, binds to ComR and activates the comX promoter. Additionally, intracellular signaling occurs with ComS binding directly to ComR. The reporter strain harbors a plasmid, pDL278, carrying a copy of gfp that is driven by the comX promoter (PcomX) to monitor ComRS signaling activation. b Cocultures of the S. mutans PcomX::gfp reporter strain grown with either S. mutans UA159 (control, green circles), S. gordonii DL1 (blue squares), S. sanguinis SK150 (orange triangles), or S. sp. A12 (red diamonds). Colored, non-connected symbols represent relative fluorescent units (RFUs) plotted on the left y-axis, while black, connected lines with symbols represent growth of the cocultures over the course of the experiment measured by optical density at 600 nm plotted on the right y-axis. Data are averages from three biological replicates of the experiment. c Percentage of each species remaining within the coculture after 18 h of monitoring, determined by colony forming unit (CFU) plating. The PcomX::gfp reporter strain is represented in the orange bars, while the competitor, listed on the left y-axis, is represented in blue. Average of collected CFUs is shown to the right. Data represent averages from three biological replicates of the experiment that was conducted in panel (b). d Cocultures of the S. mutans PcomX::gfp reporter strain in which 5 µM sXIP was added prior to the start of the experiment. e Cocultures of the S. mutans PcomX::gfp reporter strain that contains a plasmid that overexpresses the XIP peptide precursor, ComS. Control represents the PcomX::gfp reporter that contained an empty vector only.

Full size image

When the PcomX::gfp reporter strain was cocultured with wild-type S. mutans UA159 (control), robust ComRS signaling was observed as cell density increased (Fig. 1b). However, when cocultured with a competitor Streptococcus spp., no signal from the S. mutans reporter could be detected above background levels; i.e., the nonspecific fluorescence generated by an S. mutans strain that did not contain a copy of the gfp gene. The lack of fluorescence in the cocultures with commensals was not due to growth inhibition of S. mutans as the reporter strain constituted 10 ± 3%, 37 ± 5%, or 54 ± 3% of the total colony forming units (CFUs) recovered after 18 h of coculturing with S. gordonii DL1, S. sanguinis SK150, or S. sp. A12, respectively (Fig. 1c). The quantity of S. mutans cells in the commensal cocultures compared favorably with the recovery of the reporter strain (54 ± 5%) in coculture with wild-type S. mutans UA159. Of note, the fact that equal proportions of reporter and wild-type S. mutans were recovered from cocultures demonstrated that the presence of the GFP gene fusion did not compromise the fitness of the reporter strain, further verified by growth rate comparisons between wild-type and reporter strains (Supplementary Fig. 1).
Two strategies were implemented to try to recover active ComRS signaling by the reporter strain during cocultivation with commensal streptococci. First, synthetic XIP was added to the cocultures to a final concentration of 5 µM just prior to the beginning of the fluorescence monitoring phase of the experiments, and cocultures were observed as above. No detectable fluorescence signal was recorded above background in the cocultures, with or without exogenously added XIP (Fig. 1d). Second, a plasmid encoding a copy of the XIP precursor comS under the control of a highly expressed constitutive promoter (P23) [42] was introduced into the S. mutans reporter strain; we previously reported that overexpression of comS could strongly activate PcomX [28]. However, no increase in GFP expression was observed in cocultures of the comS overexpressing strain with the commensals, whereas signaling was greatly enhanced when cocultured with strain UA159 as a control (Fig. 1e).
To ensure these observations were not limited to only planktonic growth conditions, we examined S. mutans ComRS signaling in cocultured biofilm populations. While almost all cells harboring the PcomX::gfp reporter were GFP-positive in the control biofilms (coculture of the reporter with wild-type S. mutans), confocal imaging of biofilms containing competitor streptococci uniformly showed that almost no S. mutans cells were expressing detectable GFP (Fig. 2a). However, in some frames (0.22 × 0.22 mm frames, ~30,000 S. mutans cells per frame), a small number of cells (1–3 cells per frame) were GFP-positive. When 3D renderings of these areas within the biofilm were constructed, GFP-positive cells were found close to the substratum (Fig. 2b and Movie S1, same area of biofilm as top panel of Fig. 2b). Also, PcomX-active cells were not necessarily confined to distinct S. mutans microcolonies, and in some cases could be seen adjacent to the competitor streptococci, which carried a constitutively expressed red fluorescent protein (DsRed2) for their identification. To quantify the different types of cells in the biofilm populations, we physically dispersed the biofilms by sonication and analyzed the populations by flow cytometry (Supplementary Fig. 2). About 1 in 10,000 S. mutans cells counted displayed activation of PcomX within the biofilms, which was similar to the proportions of GFP-expressing cells in planktonic growth conditions (Supplemental Table 1).
Fig. 2: S. mutans peptide signaling in coculture biofilms.

a 3D volume projections of imaged biofilms in the XY-orientation (from the top looking down). Each biofilm contains either S. mutans UA159 with a constitutive gfp reporter plasmid (top row), or the PcomX::gfp reporter plasmid (bottom row) that was cocultured with either S. mutans (control; left), S. gordonii DL1 (middle), or S. sp. A12 (right) who all constitutively produce DsRed2. To the right of each expanded color image is the black and white image capture of each individual channel: blue (top), green (middle), and red (bottom). b Zoomed image frames of PcomX-active cells within cocultured biofilms with S. gordonii DL1. The images captured are a single z plane near or at the biofilm substratum. Two different areas of the biofilm (top and bottom rows) were imaged. Each panel represents one color channel of blue (SYTO 42 stained; total cells), green (PcomX::gfp positive cells), or red (S. gordonii P23::DsRed2) followed by the merged image on the far right. The top panel of (b) is the same area of biofilm shown in Movie S1.

Full size image

Commensal signaling inhibition is dependent on cell contact
Changes in phenotypes that are observed when two different species of bacteria are cocultured can usually be induced by secreted molecules from one of the bacterial strains [1]. We suspected that molecule(s) secreted by the competitor strains are required for shutting down cell–cell signaling in S. mutans. To explore this hypothesis, we cultured the competitors individually overnight and collected the supernatant fluids after centrifugation. The supernates were then filter sterilized, pH adjusted from ~6.3 to 7.0 with NaOH, and carbohydrate was added back to achieve a final concentration of added glucose to 20 mM. We then inoculated our reporter strain into the commensal supernates and monitored fluorescence activity (Fig. 3a). Surprisingly, ComRS signaling was readily observed in all supernates. In fact, reporter activity tended to be higher in the supernates of competitors compared to controls.
Fig. 3: Cell contact dependence in signaling inhibition.

a Growth and fluorescence of S. mutans PcomX::gfp reporter strain in spent supernatant fluids of either S. mutans UA159 (control, green circles), S. gordonii DL1 (blue squares), S. sanguinis SK150 (orange triangles), or S. sp. A12 (red diamonds). Depiction on top shows methods used to treat supernatant fluids following harvesting and prior to reporter strain inoculation. Overnight cultures of selected strains where centrifuged, spent supernates removed, filter sterilized, the pH was adjusted to 7.0 and 20 mM additional glucose was added. The PcomX::gfp reporter strain was then inoculated and monitored for 18 h in a Synergy 2 multimode plate reader. b Growth of cocultures in a transwell apparatus. The PcomX::gfp reporter strain was first inoculated in 0.1 mL of CDM medium in a 96-well microtiter plate. The transwell plate was then overlaid on top of the 96-well plate, and 0.1 mL of CDM inoculated with either S. mutans UA159 (control, green circles), S. gordonii DL1 (blue squares), S. sanguinis SK150 (orange triangles), or S. sp. A12 (red diamonds) was added to the top chamber, as shown. Cultures of the reporter strain and competitor were separated by a 0.4 µM pore size polycarbonate filter membrane. Fluorescence (RFUs) of the reporter strain was monitored for 18 h.

Full size image

In another experiment to confirm these results, we grew competitor and our reporter strains together in a transwell apparatus, so that both bacterial strains shared the same growth medium, but were physically separated by 0.4 µm pore size polycarbonate membrane that would allow passage of small molecule(s) between the chambers (Fig. 3b). Even in the transwell system, cell signaling was robust in cocultures containing competitor species. This result is consistent with data showing that the proximity of live commensal cells with S. mutans prior to signal activation is required for the signaling inhibition.
Impairment of S. mutans cell signaling by oral commensals is conserved across species
We next screened a collection of low-passage oral streptococci that had been previously genome sequenced [43] to determine whether the ability to inhibit S. mutans ComRS signaling was conserved across commensal species and to assess whether the presence or absence of certain genes might contribute to inhibition of peptide signaling. Ten different low-passage clinical isolates of S. gordonii, ten isolates of S. sanguinis, and five isolates of S. sp. A12-related organisms [19] were cocultured with our S. mutans ComRS signaling reporter. The S. sp. A12-related organisms included strains classified as A12-like (A13 and BCC21), as Streptococcus australis (G1 and G2), or as Streptococcus parasanguinis (A1). Interestingly, significant production of GFP by S. mutans was evident when cultured with one isolate of S. sanguinis (BCC64) and with three isolates that were classified as A12-related (BCC21, G1 and G2) (Fig. 4a). However, these results were most likely due to the inability of these isolates to grow well within the CDM medium during the course of the experiment (Supplementary Fig. 3). In fact, after 18 h of monitoring, these isolates comprised  0.1 after 12 h as monitored using a Bioscreen system, see Supplementary Fig. 3) inhibited PcomX activation. Thus, if a commensal strain could grow in CDM, even somewhat poorly, it could completely inhibit ComRS signaling.
Fig. 4: Conservation of ComRS signaling antagonism across oral isolates.

a Relative fluorescent units (RFUs) of the S. mutans PcomX::gfp reporter strain cocultured with clinical oral isolates of either S. gordonii, S. sanguinis or S. sp. A12-like strains. Relative fluorescent units were recorded after coculture inoculation at 1:1 ratio and 12 h of incubation at 37 °C. Results from four biological replicates of the experiment are shown. b RFUs after 12 h of incubation of the PcomX::gfp reporter harbored in various S. mutans clinical isolates. The PcomX::gfp reporter strain was cocultured with either S. mutans UA159 (control; black dots and bars) or an oral competitor streptococci (S. sp. A12, red dots and bars).

Full size image

We also tested several genomically and phenotypically diverse isolates of S. mutans [44, 45], both in coculture with our PcomX::gfp reporter in the UA159 background (Supplementary Fig. 4) and against competitor Streptococcus spp., after transformation of the S. mutans strains with the PcomX reporter plasmid (Fig. 4b). Various levels of spontaneous activation of the PcomX::gfp reporter were observed among the different S. mutans strains in monocultures in CDM, consistent with recent reports showing strain-dependent differences in S. mutans peptide signaling [46]. One isolate, Smu107 (R221), had undetectable levels of GFP in monoculture in CDM alone. All others showed activity above baseline. However, when cocultured with S. sp. A12, ComRS signaling was inhibited to an extent similar to that observed with strain wild-type UA159. Therefore, the ability to obstruct ComRS signaling is conserved among isolates of S. gordonii, S. sanguinis, and A12-related streptococci, and inhibition by commensals is similarly conserved in genomically diverse isolates of S. mutans.
Relatively small proportions of live commensal streptococci can inhibit signaling
To verify that the ability of the competitor species to grow (viability) was required for inhibition of peptide signaling, we used two different treatments of the competitor species S. sp. A12 after it was grown to mid-exponential phase: 80 °C for 0.5 h in a heating block (Fig. 5a) or treatment with 4% paraformaldehyde for 1 h at ambient temperature (Fig. 5b). After treatment, the inactivated commensal cells were washed and resuspended in fresh CDM and then mixed with the S. mutans reporter strain to begin the experiment. With heat-treated cells, some ComRS signal activity was evident, but not near the levels seen with S. mutans-only controls. However, when the paraformaldehyde-treated cells were used, the competitor did not inhibit signaling and fluorescence, with levels being similar to the S. mutans-only control. Importantly, we determined that there was a greater number of live cells, by plating and counting CFUs, for the competitor after heat treatment, compared to paraformaldehyde fixing (Supplementary Fig. 5), which likely explains the difference in effects on PcomX activation. These results support that metabolically active and growing competitors are required for S. mutans ComRS signaling obstruction.
Fig. 5: Importance of oral competitor cell density in signaling inhibition.

Cocultures of the S. mutans PcomX::gfp reporter strain with untreated or treated cells by either a 0.5 h heat inactivation at 80 °C or b 1 h suspension in 4% paraformaldehyde. Data represent averages from three biological replicates. c Dilution of an oral competitor streptococci (S. sp. A12) in coculture with the S. mutans PcomX::gfp reporter strain. Legend (top left) refers to the amount of S. sp. A12 within the coculture at the time of initial inoculation. Bottom: addition of either control (UA159; blue squares) or an oral competitor streptococci (S. gordonii DL1; orange triangles) at 4.5 h to a growing culture of the S. mutans PcomX::gfp strain when competence activation was d fully detected, e beginning to be detected, or f not yet detected. See Supplementary Fig. 7 for comparisons at 4.5 and 12 h, specifically.

Full size image

Based on the intermediate inhibitory effects seen with reduced proportions of a live competitor on our reporter strain, i.e. with heat-treated cells, we tested whether some minimal proportion of live competitor was required to exert effects on ComRS signaling. We utilized S. sp. A12 and varied the percentage of S. mutans and S. sp. A12 in the cocultures, after determining that the proportions of cells recovered after 18 h were similar to the proportions in the initial inocula (Supplementary Fig. 6). Complete inhibition of S. mutans ComRS signaling occurred when S. sp. A12 constituted ≥6.3% of the initial inoculum (Fig. 5c). At 3.1 or 1.6% of S. sp. A12, reporter activity was detectable, but at lower levels than when no S. sp. A12 was present. No difference in S. mutans reporter activity was observed when  More

Site description and sampling
This experimental site was established in July 2009 at the Kessler Atmospheric and Ecological Field Station (KAEFS) in the US Great Plains in McClain County, Oklahoma (34̊ 59ʹN, 97̊ 31ʹW)14,48. Experimental design and site description were described in detail previously25. Briefly, Ambrosia trifida, Solanum carolinense and Euphorbia dentate belonging to C3 forbs, and Tridens flavus, Sporobolus compositus and Sorghum halapense belonging to C4 grasses are dominant in the site25,48. Annual mean temperature is 16.3 °C and annual precipitation is 914 mm, based on Oklahoma Climatological Survey data from 1948 to 1999. The soil type of this site is Port–Pulaski–Keokuk complex with 51% of sand, 35% of silt and 13% of clay, which is a well-drained soil that is formed in loamy sediment on flood plains. The soil has a high available water holding capacity (37%), neutral pH and 1.2 g cm−3 bulk density with 1.9% total organic matter and 0.1% total nitrogen (N)25,48. Four blocks were used in the field site experiment, in which warming is a primary factor. Two levels of warming (ambient and +3 °C) were set for four pairs of 2.5 m × 1.75 m plots by utilizing a real or dummy infrared radiator (Kalglo Electronics, Bethlehem, PA, USA). In the warmed plots, a real infrared radiator was suspended 1.5 m above the ground, and the dummy infrared radiator was suspended to simulate a shading effect of the device in the control plots.
In this study, eight surface (0–15 cm) soil samples, four from the warmed and four from the control plots, were collected annually at approximately the date of peak plant biomass (September or October) from 2010 to 2016. Three soil cores (2.5 cm diameter × 15 cm depth) were taken by using a soil sampler tube in each plot and composited to have enough samples for soil chemistry, microbiology and molecular biology analyses. A total of 56 soil samples were analyzed in this study.
Environmental and soil chemical measurements
Precipitation data were obtained from the Oklahoma Mesonet Station (Washington Station)48 located 200 m away from our experiment site, and 12-month version of the standardized precipitation-evapotranspiration index (SPEI-12) was used as annual drought index49. Air temperature, soil temperature and volumetric soil water content were described in detail previously25. Specifically, air temperature and soil temperature at the depth of 7.5 cm in the center of each field plot were measured by using Constantan-copper thermocouples wired to a Campbell Scientific CR10x data logger (Campbell Scientific, UT, USA). A portable time domain reflectometer (Soil Moisture Equipment Corp.) was used to measure soil moisture from the soil surface to a 15-cm depth once or twice a month. Three measurements of soil moisture were performed in each plot and the average of three technical replicates were used in further analyses.
All soil samples were analyzed to determine soil total organic carbon (TOC), total nitrogen (TN), soil nitrate (NO3−) and ammonia (NH4+) by the Soil, Water, and Forage Analytical Laboratory at Oklahoma State University (Stillwater, OK, USA). Soil pH was measured using a pH meter with a calibrated combined glass electrode50.
Aboveground plant communities
Aboveground plant community investigations were annually conducted at peak biomass (usually September)48,51. Aboveground plant biomass, separated into C3 and C4 species, was indirectly estimated by a modified pin-touch method48,51. Detailed description of biomass estimation is provided by Sherry et al.52. A pin frame used in this study is 1 m long and have 10 pins 10 cm apart at 30° from vertical. Pins with a 0.75 m length were raised within the frame to count hits up to 1 m high (hits over 1 m are negligible at this site). The pin frame was placed in the center of each plot to record the contact numbers of the pins separately with C3 and C4 plants (e.g., leaves and stems). The contact numbers of C3 and C4 plants were then used to estimate plant biomass using calibration equations derived from calibration plots, which were located near the experimental plots. Biomass in the calibration plots was clipped at a height of 10 cm above the ground at approximately the date of peak plant biomass (September or October). All of the species in plant community within each plot were identified to estimate species richness. Clipped plant materials were oven-dried and then correlated with the total contact number. C3 and C4 plant biomasses were estimated by using the calibration equation of contact number and plant biomass. All of the species within each plot were identified to estimate species richness of plants.
Ecosystem C fluxes and soil respiration
Ecosystem C fluxes and soil respirations were measured once or twice a month between 10:00 and 15:00 (local time) from January 2010 to December 2016 by following previous methods14,48. One square aluminum frame (0.5 m × 0.5 m) was inserted in the soil at 2 cm depth in each plot to provide a flat base between the soil surface and the CO2 sampling chamber. NEE and ecosystem respiration (ER) were measured using LI-6400 portable photosynthesis system (LI-COR). Gross primary productivity (GPP) was estimated as the difference between NEE and ER. Meanwhile, soil surface respiration was monthly measured using a LI-8100A soil flux system attached to a soil CO2 flux chamber (LI-COR). Measurements were taken above a PVC collar (80 cm2 in area and 5 cm in depth) and a PVC tube (80 cm2 in area and 70 cm in depth) in each plot. The PVC tube was permanently fixed on the ground to cut off old plant roots and prevent new roots from growing inside the tube. Any aboveground parts of living plants were removed from the PVC tubes and collars before each measurement. The CO2 efflux measured above the PVC tubes represented heterotrophic respiration (Rh) from soil microbes, while that measured above the PVC collars represented soil total respiration (Rt) including heterotrophic and autotrophic respiration (Rh and Ra) from soil microbes and plant root, respectively.
Soil decomposition rate
Weighted cellulose filter paper (Whatman CAT No. 1442-090) was placed into fiberglass mesh bags and placed vertically at 0–10 cm soil depth in each plot in March 2016. All of decomposition bags were collected back in September 2016, rinsed and dried at 60 °C for weighing. The percentage of mass loss was calculated to represent soil decomposition rate.
Molecular analyses of soil samples
The C substrate utilization patterns of soil microbial communities in 2016 were analyzed by BIOLOG EcoPlateTM (BIOLOG). The BIOLOG EcoPlateTM contains 31 of the most useful labile carbon sources for soil community analysis, which are repeated three times in each plate. In this study, the plates with diluted soil supernatant (0.5 g soil with 45 mL 0.85% NaCl) were incubated in a BIOLOG OmniLog PM System at 25 °C for 4.5 days. The color change of each well was shown as absorbance curve. The net area under the absorbance versus time curve was calculated to represent physiological activity of various C sources53. The average value from three replicates was used for analyses in this study.
Soil total DNA was extracted from 1.5 g soil by freeze-grinding and SDS-based lysis54, and purified with a MoBio PowerSoil DNA isolation kit (MoBio Laboratories)25. Then, 10 ng DNA per sample were used for library construction and amplicon sequencing. Amplicons sequencing was performed with cautions in terms of experimental preparations and data analyses to ensure sequence representativeness and semi-quantitative nature55. The V4 region of bacterial and archaeal 16S rRNA genes were amplified with the primer set 515F (5ʹ-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5ʹ-GGACTACHVGGGTWTCTAAT-3ʹ), and fungal ITSs between 5.8S and 28S rRNA genes were amplified with the primer set ITS7F (5ʹ-GTGARTCATCGARTCTTTG-3ʹ) and ITS4R (5ʹ-TCCTCCGCTTATTGATATGC-3ʹ). PCR products from different samples were sequenced on a MiSeq platform (Illumina, Inc.) using 2 × 250 pair-end sequencing kit. Raw sequences were submitted to our Galaxy sequence analysis pipeline (http://zhoulab5.rccc.ou.edu:8080) to further analyze according to the protocol in the pipeline25. Finally, OTUs were clustered by UPARSE56 at 97% identity for both 16S rRNA gene and ITS. All sequences were randomly resampled to 30,000 sequences for 16S rRNA gene and 10,000 sequences for ITS per sample. Representative sequences of OTUs were annotated taxonomically by the Ribosomal Database Project (RDP) Classifier with 50% confidence estimates.
GeoChip 5.0 M, a functional gene array57, was used for all 56 samples from 2010 to 2016. GeoChip hybridization, scanning and data processing were performed in the Institute for Environmental Genomics, University of Oklahoma57,58. Specifically, 800 ng of purified soil DNA of each sample was mixed with 5.5 µl random primers (Life Technologies, random hexamers, 3 µg/µl), diluted with nuclease-free water to 35 µl, heated to 99 °C for 5 min, and placed on ice immediately. The labeling master mix (15 µl), including 0.5 µl of Cy-3 dUTP (25 nM; GE Healthcare), 2.5 µl of dNTP (2.5 mM dTTP, 5 mM dAGC-TP), 1 µl of Klenow (imer; San Diego, CA; 40 U ml−1), 5 µl Klenow buffer, and 2.5 µl of water, was added in the sample mixed solution. The samples were incubated at 37 °C for 6 h in a thermocycler, and then incubated at 95 °C for 3 min to inactivate the enzyme. Subsequently, samples were protected from the light as much as possible. Labeled DNA was cleaned using a QIAquick purification kit (Qiagen) according the manufacturer’s instructions and then dried thoroughly in a SpeedVac (45 °C, 45 min; ThermoSavant).
Labeled DNA was resuspended into 27.5 µl of DNase-free water, and then mixed completely with 99.4 µl of hybridization solution, containing 63.5 µl of formamide (10% final concentration), 2 × HI-RPM hybridization buffer, 12.7 µl of 10 × aCGH blocking agent, 0.05 μg/µl Cot-1 DNA, and 10 pM CORS58. The mixed solution was denatured at 95 °C for 3 min, and then incubated at 37 °C for 30 min. The DNA solution was centrifuged at 6000 × g for 1 min to collect liquid at the bottom of the tube. 110 µl of the solution was pipetted into the center of the well of the gasket slide. The array slide was placed on the gasket slide, sealed using a SureHyb chamber, hybridized at 67 °C for 24 h at 20 rpm in a hybridization oven. After hybridization, slides were washed in room temperature with Wash Buffer 1 (Agilent) and Wash Buffer 2 (Agilent).
The slides were imaged as a Multi-TIFF with a NimbleGen MS200 Microarray Scanner (Roche NimbleGen, Inc., Madison, WI, United States). The raw signals from NimbleGen were submitted to the Microarray Data Manager on our website (http://ieg.ou.edu/microarray), cleaned, normalized and analyzed using the data-analysis pipeline. Briefly, probe quality was assessed, and poor or low signal probes were removed. Probe spots with coefficient of variance (CV; probe signal SD/signal) >0.8 were removed. Then, the signal-to-noise ratio (SNR) was calculated. As suggested by Agilent, the average signal of Agilent’s negative control probes within each subarray was used as the background signal for the probes in that subarray instead of the local background typically used. The signal intensity for each spot was corrected by subtracting the background signal intensity. If the net difference was30%, aligned length >20 a.a., and e-value More