Philippa Kaur

1.
Faridah-Hanum, I. et al. Development of a comprehensive mangrove quality index (MQI) in Matang Mangrove: assessing mangrove ecosystem health. Ecol. Ind. 102, 103–117 (2019).
Article  Google Scholar 
2.
Spalding, M. World Atlas of Mangroves (Routledge, New York, 2010).
Google Scholar 

3.
Himes-Cornell, A., Grose, S. O. & Pendleton, L. Mangrove ecosystem service values and methodological approaches to valuation: where do we stand?. Front. Mar. Sci. 5, 376 (2018).
Article  Google Scholar 

4.
Wu, Y., Ricklefs, R. E., Huang, Z., Zan, Q. & Yu, S. Winter temperature structures mangrove species distributions and assemblage composition in China. Glob. Ecol. Biogeogr. 27, 1492–1506 (2018).
Article  Google Scholar 

5.
Das, S. S., Das, S. & Ghosh, P. Phylogenetic relationships among the mangrove species of Acanthaceae found in Indian Sundarban, as revealed by RAPD analysis. Adv. Appl. Sci. Res 6, 179–184 (2015).
CAS  Google Scholar 

6.
Kathiresan, K. & Bingham, B. L. Biology of mangroves and mangrove ecosystems. Adv. Mar. Biol. 40, 84–254 (2001).
Google Scholar 

7.
Sannigrahi, S. et al. Responses of ecosystem services to natural and anthropogenic forcings: a spatial regression based assessment in the world’s largest mangrove ecosystem. Sci. Total Environ. 715, 137004 (2020).
ADS  CAS  PubMed  Article  Google Scholar 

8.
dos Santos, N. M. & Lana, P. Present and past uses of mangrove wood in the subtropical Bay of Paranaguá (Paraná, Brazil). Ocean Coast. Manag. 148, 97–103 (2017).
Article  Google Scholar 

9.
Jusoff, K. Malaysian mangrove forests and their significance to the coastal marine environment. Pol. J. Environ. Stud. 22, 979–1005 (2013).
Google Scholar 

10.
Duke, N. C., Lo, E. & Sun, M. Global distribution and genetic discontinuities of mangroves–emerging patterns in the evolution of Rhizophora. Trees 16, 65–79 (2002).
Article  Google Scholar 

11.
Nagelkerken, I. et al. The habitat function of mangroves for terrestrial and marine fauna: a review. Aquat. Bot. 89, 155–185 (2008).
Article  Google Scholar 

12.
Tomlinson, P. B. The Botany of Mangroves (Cambridge University Press, Cambridge, 2016).
Google Scholar 

13.
Dahdouh-Guebas, F. et al. How effective were mangroves as a defence against the recent tsunami?. Curr. Biol. 15, R443–R447 (2005).
CAS  PubMed  Article  PubMed Central  Google Scholar 

14.
Saenger, P. Mangrove Ecology, Silviculture and Conservation (Springer, Berlin, 2002).
Google Scholar 

15.
Lakshmi, M., Parani, M. & Parida, A. Molecular phylogeny of mangroves IX: molecular marker assisted intra-specific variation and species relationships in the Indian mangrove tribe Rhizophoreae. Aquat. Bot. 74, 201–217 (2002).
CAS  Article  Google Scholar 

16.
Grover, A. & Sharma, P. Development and use of molecular markers: past and present. Crit. Rev. Biotechnol. 36, 290–302 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

17.
Adsul, G. G., Chaurasia, A. K., Dhake, A. V. & Kothari, R. M. RAPD analysis of phylogenetic relationships and genetic variations in genus Allium. Biochem. Indian J. 3, 1–5 (2015).
Google Scholar 

18.
Xu, K. et al. Identification of tuna species (Thunnini tribe) by PCR-RFLP analysis of mitochondrial DNA fragments. Food Agric. Immunol. 27, 301–313 (2016).
CAS  Article  Google Scholar 

19.
Wambugu, P. W., Brozynska, M., Furtado, A., Waters, D. L. & Henry, R. J. Relationships of wild and domesticated rices (Oryza AA genome species) based upon whole chloroplast genome sequences. Sci. Rep. 5, 1–9 (2015).
Article  Google Scholar 

20.
Middleton, C. P. et al. Sequencing of chloroplast genomes from wheat, barley, rye and their relatives provides a detailed insight into the evolution of the Triticeae tribe. PLoS ONE 9, e85761 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

21.
Raman, G. & Park, S. The complete chloroplast genome sequence of Ampelopsis: gene organization, comparative analysis, and phylogenetic relationships to other angiosperms. Front. Plant Sci. 7, 341 (2016).
PubMed  PubMed Central  Article  Google Scholar 

22.
Yang, J. B., Li, D. Z. & Li, H. T. Highly effective sequencing whole chloroplast genomes of angiosperms by nine novel universal primer pairs. Mol. Ecol. Resour. 14, 1024–1031 (2014).
CAS  PubMed  Article  PubMed Central  Google Scholar 

23.
Su, H.-J., Hogenhout, S. A., Al-Sadi, A. M. & Kuo, C.-H. Complete chloroplast genome sequence of Omani lime (Citrus aurantiifolia) and comparative analysis within the rosids. PLoS ONE 9, e113049 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

24.
Hu, S. Phylogeny and Chloroplast Evolution in BRASSICACEAE, University of Trento, (2016).

25.
Santos, C. G. Development of new tools for the identification of plants using chloroplast DNA sequences. (2018).

26.
Singh, B. P., Kumar, A., Kaur, H., Singh, H. & Nagpal, A. K. CpGDB: A comprehensive database of chloroplast genomes. Bioinformation 16, 171 (2020).
PubMed  PubMed Central  Article  Google Scholar 

27.
Wu, Z. The new completed genome of purple willow (Salix purpurea) and conserved chloroplast genome structure of Salicaceae. J. Nat. Sci 1, e49 (2015).
Google Scholar 

28.
Wu, Z. The whole chloroplast genome of shrub willows (Salix suchowensis). Mitochondrial DNA Part A 27, 2153–2154 (2016).
CAS  Google Scholar 

29.
Egamberdiev, S. S. et al. Comparative assessment of genetic diversity in cytoplasmic and nuclear genome of upland cotton. Genetica 144, 289–306 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

30.
Williams, A. V., Miller, J. T., Small, I., Nevill, P. G. & Boykin, L. M. Integration of complete chloroplast genome sequences with small amplicon datasets improves phylogenetic resolution in Acacia. Mol. Phylogenet. Evol. 96, 1–8 (2016).
CAS  PubMed  Article  Google Scholar 

31.
Fučíková, K. et al. New phylogenetic hypotheses for the core Chlorophyta based on chloroplast sequence data. Front. Ecol. Evol. 2, 63 (2014).
Google Scholar 

32.
Friis, G. et al. A high-quality genome assembly and annotation of the gray mangrove, Avicennia marina. bioRxiv (2020).

33.
Asaf, S. et al. Complete chloroplast genome of Nicotiana otophora and its comparison with related species. Front. Plant Sci. 7, 843 (2016).
PubMed  PubMed Central  Article  Google Scholar 

34.
Asaf, S., Khan, A. L., Khan, A. & Al-Harrasi, A. Unraveling the chloroplast genomes of two prosopis species to identify its genomic information, comparative analyses and phylogenetic relationship. Int. J. Mol. Sci. 21, 3280 (2020).
CAS  PubMed Central  Article  PubMed  Google Scholar 

35.
Souladeth, P., Tagane, S., Zhang, M., Okabe, N. & Yahara, T. Flora of Nam Kading National Protected Area I: a new species of yellow-flowered Strobilanthes (Acanthaceae), S. namkadingensis. PhytoKeys 81, 11 (2017).
Article  Google Scholar 

36.
Yang, L. et al. The complete chloroplast genome of Swertia tetraptera and phylogenetic analysis. Mitochondrial DNA Part B 5, 164–165 (2020).
Article  Google Scholar 

37.
Biju, V. C. et al. The complete chloroplast genome of Trichopus zeylanicus, and phylogenetic analysis with dioscoreales. Plant Genome 12, 190032 (2019).
CAS  Article  Google Scholar 

38.
Silva, S. R. et al. The chloroplast genome of Utricularia reniformis sheds light on the evolution of the ndh gene complex of terrestrial carnivorous plants from the Lentibulariaceae family. PLoS ONE 11, e0165176 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

39.
Zuo, L.-H. et al. The first complete chloroplast genome sequences of Ulmus species by de novo sequencing: Genome comparative and taxonomic position analysis. PLoS ONE 12, e0171264 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

40.
Niu, Z., Huang, S., Deng, Y. & Chen, X. The complete chloroplast genome of Justicia leptostachya (Acanthaceae). Mitochondrial DNA Part B 4, 1114–1115 (2019).
Article  Google Scholar 

41.
Chen, H. et al. Sequencing and analysis of Strobilanthes cusia (Nees) Kuntze chloroplast Genome revealed the rare simultaneous contraction and expansion of the inverted repeat region in Angiosperm. Front. Plant Sci. 9, 324 (2018).
PubMed  PubMed Central  Article  Google Scholar 

42.
Ding, P. et al. The complete chloroplast genome sequence of the medicinal plant Andrographis paniculata. Mitochondrial DNA Part A 27, 2347–2348 (2016).
CAS  Article  Google Scholar 

43.
Yaradua, S. S., Alzahrani, D. A., Albokhary, E. J., Abba, A. & Bello, A. Complete chloroplast genome sequence of Justicia flava: genome comparative analysis and phylogenetic relationships among Acanthaceae. BioMed Res. Int. 2019, 4370258 (2019).
PubMed  PubMed Central  Article  CAS  Google Scholar 

44.
Huang, S., Deng, Y. & Ge, X. The complete chloroplast genome of Aphelandra knappiae (Acanthaceae). Mitochondrial DNA Part B 4, 273–274 (2019).
Article  Google Scholar 

45.
Li, M.-N. et al. Complete plastome sequence of Clinacanthus nutans (Acanthaceae): a medicinal species in Southern China. Mitochondrial DNA Part B 4, 118–119 (2019).
Article  Google Scholar 

46.
Jiang, M., Wang, J. & Zhang, H. Characterization and phylogenetic analysis of the complete chloroplast genome sequence of Disanthus cercidifolius subsp. longipes (Hamamelidaceae), a rare and endangered wild plant species in China. Mitochondrial DNA Part B 5, 1206–1207 (2020).
CAS  PubMed  Article  Google Scholar 

47.
Steane, D. A. Complete nucleotide sequence of the chloroplast genome from the Tasmanian blue gum, Eucalyptus globulus (Myrtaceae). DNA Res. 12, 215–220 (2005).
CAS  PubMed  Article  Google Scholar 

48.
Park, J., Kim, Y., Xi, H. & Heo, K.-I. The complete chloroplast genome of ornamental coffee tree, Coffea arabica L. (Rubiaceae). Mitochondrial DNA Part B 4, 1059–1060 (2019).
Article  Google Scholar 

49.
Park, J. et al. The complete chloroplast genome of common camellia tree, Camellia japonica L. (Theaceae), adapted to cold environment in Korea. Mitochondrial DNA Part B 4, 1038–1040 (2019).
Article  Google Scholar 

50.
Arif Khan, S. A. et al. First complete chloroplast genomics and comparative phylogenetic analysis of Commiphora gileadensis and C. foliacea: Myrrh producing trees. PLoS ONE 14, e0225469 (2019).
PubMed  PubMed Central  Article  CAS  Google Scholar 

51.
Xu, J. et al. The first intron of rice EPSP synthase enhances expression of foreign gene. Sci. China Ser. C Life Sci. 46, 561 (2003).
CAS  Article  Google Scholar 

52.
Kelchner, S. A. The evolution of non-coding chloroplast DNA and its application in plant systematics. Ann. Missouri Bot. Gard. 87, 482–498 (2000).
Article  Google Scholar 

53.
Khan, A. et al. Complete chloroplast genomes of medicinally important Teucrium species and comparative analyses with related species from Lamiaceae. PeerJ 7, e7260 (2019).
PubMed  PubMed Central  Article  Google Scholar 

54.
Asaf, S., Khan, A., Khan, A. L., Al-Harrasi, A. & Al-Rawahi, A. Complete chloroplast genomes of Vachellia nilotica and Senegalia senegal: comparative genomics and phylogenomic placement in a new generic system. PLoS ONE 14, e0225469 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

55.
Su, Y., He, Z., Wang, Z., Hong, Y. & Wang, T. Characterization of the complete chloroplast genome of Leptochilus decurrens (Polypodiaceae), a least concern folk medicinal fern. Mitochondrial DNA Part B 4, 3346–3347 (2019).
PubMed  Article  PubMed Central  Google Scholar 

56.
Kumar, S. & Shanker, A. In silico comparative analysis of simple sequence repeats in chloroplast genomes of genus nymphaea. J. Sci. Res. 64, 186–192 (2020).
Google Scholar 

57.
Asaf, S. et al. The complete chloroplast genome of wild rice (Oryza minuta) and its comparison to related species. Front. Plant Sci. 8, 304 (2017).
PubMed  PubMed Central  Article  Google Scholar 

58.
Asaf, S. et al. Complete chloroplast genome sequence and comparative analysis of loblolly pine (Pinus taeda L.) with related species. PLoS ONE 13, e0192966 (2018).
PubMed  PubMed Central  Article  CAS  Google Scholar 

59.
Khan, A. L. et al. The first chloroplast genome sequence of Boswellia sacra, a resin-producing plant in Oman. PLoS ONE 12, e0169794 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

60.
Asaf, S. et al. Comparative analysis of complete plastid genomes from wild soybean (Glycine soja) and nine other Glycine species. PLoS ONE 12, e0182281 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

61.
Qian, J. et al. The complete chloroplast genome sequence of the medicinal plant Salvia miltiorrhiza. PLoS ONE 8, e57607 (2013).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

62.
Wang, W. & Messing, J. High-throughput sequencing of three Lemnoideae (duckweeds) chloroplast genomes from total DNA. PLoS ONE 6, e24670 (2011).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

63.
Shen, X. et al. Complete chloroplast genome sequence and phylogenetic analysis of the medicinal plant Artemisia annua. Molecules 22, 1330 (2017).
ADS  PubMed Central  Article  CAS  Google Scholar 

64.
Liu, M. et al. The complete chloroplast genome sequence of Tartary Buckwheat Cultivar Miqiao 1 (Fagopyrum tataricum Gaertn.). Mitochondrial DNA Part B 1, 577–578 (2016).
PubMed  Article  PubMed Central  Google Scholar 

65.
Fu, P.-C., Zhang, Y.-Z., Geng, H.-M. & Chen, S.-L. The complete chloroplast genome sequence of Gentiana lawrencei var. farreri (Gentianaceae) and comparative analysis with its congeneric species. PeerJ 4, e2540 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

66.
Choi, K. S., Chung, M. G. & Park, S. The complete chloroplast genome sequences of three Veroniceae species (Plantaginaceae): comparative analysis and highly divergent regions. Front. Plant Sci. 7, 355 (2016).
PubMed  PubMed Central  Google Scholar 

67.
Cheon, K.-S., Kim, K.-A., Kwak, M., Lee, B. & Yoo, K.-O. The complete chloroplast genome sequences of four Viola species (Violaceae) and comparative analyses with its congeneric species. PLoS ONE 14, e0214162 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

68.
Munyao, J. N. et al. Complete chloroplast genomes of chlorophytum comosum and chlorophytum gallabatense: genome structures, comparative and phylogenetic analysis. Plants 9, 296 (2020).
CAS  PubMed Central  Article  Google Scholar 

69.
Yu, T., Huang, B.-H., Zhang, Y., Liao, P.-C. & Li, J.-Q. Chloroplast genome of an extremely endangered conifer Thuja sutchuenensis Franch.: gene organization, comparative and phylogenetic analysis. Physiol. Mol. Biol. Plants 1–10 (2020).

70.
Sabri, D. M., El-Hussieny, S. A. & Elnwishy, N. Genotypic Variations of Mangrove (Avicennia marina) in Nabq Protectorate, South Sinai Egypt. Int. J. Agric. Biol. 20, 637–646 (2018).
Article  Google Scholar 

71.
71Basyuni, M., Baba, S. & Oku, H. in IOP Conference Series: Materials Science and Engineering. 012247 (IOP Publishing).

72.
Sahu, S. K., Singh, R. & Kathiresan, K. Multi-gene phylogenetic analysis reveals the multiple origin and evolution of mangrove physiological traits through exaptation. Estuar. Coast. Shelf Sci. 183, 41–51 (2016).
ADS  Article  Google Scholar 

73.
Khan, A. et al. First complete chloroplast genomics and comparative phylogenetic analysis of Commiphora gileadensis and C. foliacea: Myrrh producing trees. PLoS ONE 14, e0208511 (2019).
CAS  PubMed  PubMed Central  Article  Google Scholar 

74.
Arif, M. F., Aristya, G. R., Subositi, D., Sari, A. N. & Kasiamdari, R. S. rbcL and matK chloroplast DNA composition of green chireta (Andrographis paniculata) from Indonesia. Biodivers. J. Biol. Divers. 20, 3575–3583 (2019).
Google Scholar 

75.
Langdon, W. B. Performance of genetic programming optimised Bowtie2 on genome comparison and analytic testing (GCAT) benchmarks. BioData Min. 8, 1 (2015).
MathSciNet  CAS  PubMed  PubMed Central  Article  Google Scholar 

76.
Kearse, M. et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28, 1647–1649 (2012).
PubMed  PubMed Central  Article  Google Scholar 

77.
Gan, H. M., Schultz, M. B. & Austin, C. M. Integrated shotgun sequencing and bioinformatics pipeline allows ultra-fast mitogenome recovery and confirms substantial gene rearrangements in Australian freshwater crayfishes. BMC Evol. Biol. 14, 19 (2014).
PubMed  PubMed Central  Article  CAS  Google Scholar 

78.
Feng, G., Yang, J. & Peng, F.-R. Characterization of complete chloroplast genome of artificial hybrid passion fruit ‘Ziyan’, Passiflora edulis Sims× P. edulis f. edulis Sims (Passifloraceae). Mitochondrial DNA Part B 5, 1720–1721 (2020).
Article  Google Scholar 

79.
Brown, J., Pirrung, M. & McCue, L. A. FQC Dashboard: integrates FastQC results into a web-based, interactive, and extensible FASTQ quality control tool. Bioinformatics 33, 3137–3139 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 

80.
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).
CAS  PubMed  PubMed Central  Article  Google Scholar 

81.
Wyman, S. K., Jansen, R. K. & Boore, J. L. Automatic annotation of organellar genomes with DOGMA. Bioinformatics 20, 3252–3255 (2004).
CAS  PubMed  PubMed Central  Article  Google Scholar 

82.
Schattner, P., Brooks, A. N. & Lowe, T. M. The tRNAscan-SE, snoscan and snoGPS web servers for the detection of tRNAs and snoRNAs. Nucleic Acids Res. 33, W686–W689 (2005).
CAS  PubMed  PubMed Central  Article  Google Scholar 

83.
Lohse, M., Drechsel, O. & Bock, R. OrganellarGenomeDRAW (OGDRAW): a tool for the easy generation of high-quality custom graphical maps of plastid and mitochondrial genomes. Curr. Genet. 52, 267–274 (2007).
CAS  PubMed  Article  PubMed Central  Google Scholar 

84.
Kumar, S., Stecher, G., Li, M., Knyaz, C. & Tamura, K. MEGA X: molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35, 1547–1549 (2018).
CAS  PubMed  PubMed Central  Article  Google Scholar 

85.
Frazer, K. A., Pachter, L., Poliakov, A., Rubin, E. M. & Dubchak, I. VISTA: computational tools for comparative genomics. Nucleic Acids Res. 32, W273–W279 (2004).
CAS  PubMed  PubMed Central  Article  Google Scholar 

86.
Kurtz, S. et al. REPuter: the manifold applications of repeat analysis on a genomic scale. Nucleic Acids Res. 29, 4633–4642 (2001).
MathSciNet  CAS  PubMed  PubMed Central  Article  Google Scholar 

87.
Beier, S., Thiel, T., Münch, T., Scholz, U. & Mascher, M. MISA-web: a web server for microsatellite prediction. Bioinformatics 33, 2583–2585 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 

88.
Wirawan, A., Kwoh, C. K., Hsu, L. Y. & Koh, T. H. in International Conference on Computational Systems-Biology and Bioinformatics. 151–164 (Springer).

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

90.
Srivathsan, A. & Meier, R. On the inappropriate use of Kimura-2-parameter (K2P) divergences in the DNA-barcoding literature. Cladistics 28, 190–194 (2012).
Article  Google Scholar 

91.
Kumar, S., Nei, M., Dudley, J. & Tamura, K. MEGA: a biologist-centric software for evolutionary analysis of DNA and protein sequences. Brief. Bioinf. 9, 299–306 (2008).
CAS  Article  Google Scholar 

92.
Swofford, D. L. Paup*: Phylogenetic analysis using parsimony (and other methods) 4.0. B5. (2001).

93.
Wu, Z., Tembrock, L. R. & Ge, S. Are differences in genomic data sets due to true biological variants or errors in genome assembly: an example from two chloroplast genomes. PLoS ONE 10, e0118019 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar  More

1.
Ceballos, G. et al. Accelerated modern human-induced species losses: entering the sixth mass extinction. Sci. Adv. 1, e1400253 (2015).
PubMed  PubMed Central  Article  Google Scholar 
2.
De Groot, R. S., Alkemade, R., Braat, L., Hein, L. & Willemen, L. Challenges in integrating the concept of ecosystem services and values in landscape planning, management and decision making. Ecol. Complex. 7, 260–272 (2010).
Article  Google Scholar 

3.
Tilman, D. et al. Future threats to biodiversity and pathways to their prevention. Nature 546, 73–81 (2017).
CAS  PubMed  Article  PubMed Central  Google Scholar 

4.
Protected Planet Report 2016 (UNEP-WCMC and IUCN, 2016).

5.
Barnes, M. D., Glew, L., Wyborn, C. & Craigie, I. D. Prevent perverse outcomes from global protected area policy. Nat. Ecol. Evol. 2, 759–762 (2018).
PubMed  Article  PubMed Central  Google Scholar 

6.
Jones, K. R. et al. One-third of global protected land is under intense human pressure. Science 360, 788–791 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

7.
Geldmann, J. et al. Effectiveness of terrestrial protected areas in reducing habitat loss and population declines. Biol. Conserv. 161, 230–238 (2013).
Article  Google Scholar 

8.
The State of the World’s Forests 2020 (FAO and UNEP, 2020).

9.
Betts, M. G. et al. Global forest loss disproportionately erodes biodiversity in intact landscapes. Nature 547, 441–444 (2017).
CAS  PubMed  Article  PubMed Central  Google Scholar 

10.
Griscom, B. W. et al. Natural climate solutions. Proc. Natl Acad. Sci. USA 114, 11645–11650 (2017).
CAS  PubMed  Article  PubMed Central  Google Scholar 

11.
Gray, C. L. et al. Local biodiversity is higher inside than outside terrestrial protected areas worldwide. Nat. Commun. 7, 12306 (2016).
CAS  PubMed  PubMed Central  Article  Google Scholar 

12.
Coetzee, B. W., Gaston, K. J. & Chown, S. L. Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta-analysis. PLoS ONE 9, e105824 (2014).
PubMed  PubMed Central  Article  CAS  Google Scholar 

13.
Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover change. Science 342, 850–853 (2013).
CAS  Article  Google Scholar 

14.
Nelson, A. & Chomitz, K. M. Effectiveness of strict vs. multiple use protected areas in reducing tropical forest fires: a global analysis using matching methods. PLoS ONE 6, e22722 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

15.
Nolte, C., Agrawal, A., Silvius, K. M. & Soares-Filho, B. S. Governance regime and location influence avoided deforestation success of protected areas in the Brazilian Amazon. Proc. Natl Acad. Sci. USA 110, 4956–4961 (2013).
CAS  PubMed  Article  PubMed Central  Google Scholar 

16.
Spracklen, B., Kalamandeen, M., Galbraith, D., Gloor, E. & Spracklen, D. V. A global analysis of deforestation in moist tropical forest protected areas. PLoS ONE 10, e0143886 (2015).
CAS  PubMed  PubMed Central  Article  Google Scholar 

17.
Geldmann, J., Manica, A., Burgess, N. D., Coad, L. & Balmford, A. A global-level assessment of the effectiveness of protected areas at resisting anthropogenic pressures. Proc. Natl Acad. Sci. USA 116, 23209–23215 (2019).
CAS  PubMed  Article  PubMed Central  Google Scholar 

18.
Ewers, R. M. & Rodrigues, A. S. Estimates of reserve effectiveness are confounded by leakage. Trends Ecol. Evol. 23, 113–116 (2008).
PubMed  Article  PubMed Central  Google Scholar 

19.
Fuller, C., Ondei, S., Brook, B. W. & Buettel, J. C. First, do no harm: a systematic review of deforestation spillovers from protected areas. Glob. Ecol. Conserv. 18, e00591 (2019).
Article  Google Scholar 

20.
Stolton, S. et al. in Protected Area Governance and Management (eds Worboys, G. L. et al.) 145–168 (ANU Press, 2015).

21.
Scharlemann, J. P. et al. Securing tropical forest carbon: the contribution of protected areas to REDD. Oryx 44, 352–357 (2010).
Article  Google Scholar 

22.
Barnes, M. D. et al. Wildlife population trends in protected areas predicted by national socio-economic metrics and body size. Nat. Commun. 7, 12747 (2016).
CAS  PubMed  PubMed Central  Article  Google Scholar 

23.
Geldmann, J. et al. A global analysis of management capacity and ecological outcomes in terrestrial protected areas. Conserv. Lett. 11, e12434 (2018).
Article  Google Scholar 

24.
Amano, T. et al. Successful conservation of global waterbird populations depends on effective governance. Nature 553, 199–202 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

25.
Leader-Williams, N. & Albon, S. Allocation of resources for conservation. Nature 336, 533–535 (1988).
Article  Google Scholar 

26.
Jachmann, H. Monitoring law-enforcement performance in nine protected areas in Ghana. Biol. Conserv. 141, 89–99 (2008).
Article  Google Scholar 

27.
Critchlow, R. et al. Improving law-enforcement effectiveness and efficiency in protected areas using ranger-collected monitoring data. Conserv. Lett. 10, 572–580 (2017).
Article  Google Scholar 

28.
Coad, L. et al. Widespread shortfalls in protected area resourcing undermine efforts to conserve biodiversity. Front. Ecol. Environ. 17, 259–264 (2019).
Article  Google Scholar 

29.
Waldron, A. et al. Targeting global conservation funding to limit immediate biodiversity declines. Proc. Natl Acad. Sci. USA 110, 12144–12148 (2013).
CAS  PubMed  Article  PubMed Central  Google Scholar 

30.
Watson, J. E., Dudley, N., Segan, D. B. & Hockings, M. The performance and potential of protected areas. Nature 515, 67–73 (2014).
CAS  PubMed  Article  PubMed Central  Google Scholar 

31.
Bruner, A. G., Gullison, R. E. & Balmford, A. Financial costs and shortfalls of managing and expanding protected-area systems in developing countries. BioScience 54, 1119–1126 (2004).
Article  Google Scholar 

32.
Pan, Y. et al. A large and persistent carbon sink in the world’s forests. Science 333, 988–993 (2011).

33.
Report of the Conference of the Parties on its Sixteenth Session, held in Cancun from 29 November to 10 December 2010. Addendum. Part Two: Action Taken by the Conference of the Parties at its Sixteenth Session Report FCCC/CP/2010/7/Add.1 (UNFCCC, 2011).

34.
Fletcher, R., Dressler, W., Büscher, B. & Anderson, Z. R. Questioning REDD+ and the future of market-based conservation. Conserv. Biol. 30, 673–675 (2016).
PubMed  Article  Google Scholar 

35.
Ministerio de Ambiente y Desarrollo Sostenible, Instituto de Investigación de Recursos Biológicos Alexander von Humboldt Política Nacional para la Gestión Integral de la Biodiversidad y Sus Servicios Ecosistémicos (MADS, 2012).

36.
Sims, K. R. E. & Alix-Garcia, J. M. Parks versus PES: evaluating direct and incentive-based land conservation in Mexico. J. Environ. Econ. Manag. 86, 8–28 (2017).
Article  Google Scholar 

37.
James, A. N., Green, M. J. B. & Paine, J. R. A Global Review of Protected Area Budgets and Staff WCMC Biodiversity Series No.10 (World Conservation Press, 1999).

38.
Walker, S., Price, R., Rutledge, D., Stephens, R. T. & Lee, W. G. Recent loss of indigenous cover in New Zealand. New Zeal. J. Ecol. 30, 169–177 (2006).
Google Scholar 

39.
Ewers, R. M. et al. Past and future trajectories of forest loss in New Zealand. Biol. Conserv. 133, 312–325 (2006).
Article  Google Scholar 

40.
Sodhi, N. S. et al. The state and conservation of Southeast Asian biodiversity. Biodivers. Conserv. 19, 317–328 (2010).
Article  Google Scholar 

41.
Grossman, G. M. & Krueger, A. B. Environmental Impacts of a North American Free Trade Agreement (National Bureau of Economic Research, 1991).

42.
Locke, H. et al. Three global conditions for biodiversity conservation and sustainable use: an implementation framework. Natl Sci. Rev. 6, 1080–1082 (2019).
Article  Google Scholar 

43.
Lenzen, M. et al. International trade drives biodiversity threats in developing nations. Nature 486, 109–112 (2012).
CAS  PubMed  Article  PubMed Central  Google Scholar 

44.
Walker, N., Patel, S., Davies, F., Milledge, S. & Hulse, J. Demand-Side Interventions to Reduce Deforestation and Forest Degradation (International Institute for Environment and Development, 2013).

45.
Marie-Vivien, D., Garcia, C. A., Kushalappa, C. G. & Vaast, P. Trademarks, geographical indications and environmental labelling to promote biodiversity: the case of agroforestry coffee in India. Dev. Policy Rev. 32, 379–398 (2014).
Google Scholar 

46.
Symes, W. S., Rao, M., Mascia, M. B. & Carrasco, L. R. Why do we lose protected areas? Factors influencing protected area downgrading, downsizing and degazettement in the tropics and subtropics. Glob. Change Biol. 22, 656–665 (2016).
Article  Google Scholar 

47.
Adams, W. M. et al. Biodiversity conservation and the eradication of poverty. Science 306, 1146–1149 (2004).
CAS  PubMed  Article  PubMed Central  Google Scholar 

48.
Belle, E. et al. Protected Planet Report 2018 (UNEP-WCMC, IUCN and NGS, 2018).

49.
Geldmann, J. et al. Essential indicators for measuring area-based conservation effectiveness in the post-2020 global biodiversity framework. Preprint at https://doi.org/10.20944/preprints202003.0370.v1 (2020).

50.
Protected Areas Management Effectiveness Methodologies (Protected Planet 2020); http://go.nature.com/3ptIPHA

51.
Ervin, J. Rapid assessment of protected area management effectiveness in four countries. BioScience 53, 833–841 (2003).
Article  Google Scholar 

52.
Conservancy, N. Conservation Action Planning: Developing Strategies, Taking Action, and Measuring Success at any Scale: Overview of Basic Practices (Nature Conservancy, 2007).

53.
Hockings, M. et al. The World Heritage Management Effectiveness Workbook: 2007 Edition: How to Build Monitoring, Assessment and Reporting Systems to Improve the Management Effectiveness of Natural World Heritage Sites 3rd draft (Univ. Queensland, 2007).

54.
Moomaw, W. R., Masino, S. A. & Faison, E. K. Intact forests in the United States: proforestation mitigates climate change and serves the greatest good. Front. For. Glob. Change 2, 27 (2019).
Article  Google Scholar 

55.
Stolton, S., Hockings, M., Dudley, N., MacKinnon, K. & Whitten, T. Reporting Progress in Protected Areas: A Site-Level Management Effectiveness Tracking Tool (World Bank/WWF Alliance for Forest Conservation and Sustainable Use, 2003).

56.
Hockings, M. et al. The IUCN green list of protected and conserved areas: setting the standard for effective area-based conservation. Parks 25, 57–66 (2019).
Google Scholar 

57.
Locke, H. Nature needs half: a necessary and hopeful new agenda for protected areas. Nat. N. South Wales 58, 7–17 (2014).
Google Scholar 

58.
Wilson, E. O. Half-Earth: Our Planet’s Fight for Life (WW Norton & Company, 2016).

59.
The World Database on Protected Areas (WDPA) (IUCN and UNEP-WCMC, accessed 1 January 2020); https://www.protectedplanet.net/

60.
Iacus, S. M., King, G. & Porro, G. Causal inference without balance checking: coarsened exact matching. Polit. Anal. 20, 1–24 (2012).
Article  Google Scholar 

61.
Stuart, E. A. Matching methods for causal inference: a review and a look forward. Stat. Sci. 25, 1–21 (2010).
PubMed  PubMed Central  Article  Google Scholar 

62.
Schleicher, J. et al. Statistical matching for conservation science. Conserv. Biol. 34, 538–549 (2019).
PubMed  PubMed Central  Article  Google Scholar 

63.
Weiss, D. J. et al. A global map of travel time to cities to assess inequalities in accessibility in 2015. Nature 553, 333–336 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

64.
Gridded Population of the World, Version 4 (GPWv4): Population Density, Revision 11 (Columbia Univ. Center for International Earth Science Information Network, 2018).

65.
Olson, D. M. et al. Terrestrial ecoregions of the world: a new map of life on earth a new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity. BioScience 51, 933–938 (2001).
Article  Google Scholar 

66.
Curtis, P. G., Slay, C. M., Harris, N. L., Tyukavina, A. & Hansen, M. C. Classifying drivers of global forest loss. Science 361, 1108–1111 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

67.
Bode, M., Tulloch, A. I., Mills, M., Venter, O. & Ando, W. A. A conservation planning approach to mitigate the impacts of leakage from protected area networks. Conserv. Biol. 29, 765–774 (2015).
PubMed  Article  PubMed Central  Google Scholar 

68.
Carranza, T., Balmford, A., Kapos, V. & Manica, A. Protected area effectiveness in reducing conversion in a rapidly vanishing ecosystem: the Brazilian Cerrado. Conserv. Lett. 7, 216–223 (2014).
Article  Google Scholar 

69.
Ferraro, P. J. Counterfactual thinking and impact evaluation in environmental policy. New Dir. Eval. 2009, 75–84 (2009).
Article  Google Scholar 

70.
Joppa, L. N. & Pfaff, A. Global protected area impacts. Proc. R. Soc. B 278, 1633–1638 (2010).
PubMed  Article  PubMed Central  Google Scholar 

71.
Iacus, S. M., King, G. & Porro, G. CEM: software for coarsened exact matching. J. Stat. Softw. 30, 1–27 (2009).
Article  Google Scholar 

72.
Rosenbaum, P. R. Sensitivity analysis for m-estimates, tests, and confidence intervals in matched observational studies. Biometrics 63, 456–464 (2007).
CAS  PubMed  Article  PubMed Central  Google Scholar 

73.
Keele, L. An Overview of rbounds: an R Package for Rosenbaum Bounds Sensitivity Analysis with Matched Data White Paper, Columbus 1–15 (2010); https://go.nature.com/2M5DKXM

74.
Keele, L. J. rbounds: Perform Rosenbaum Bounds Sensitivity Tests for Matched and Unmatched Data. R Package (2014); https://cran.r-project.org/package=rbounds

75.
World Development Indicators 2018 (World Bank, 2018).

76.
Conner, M. M., Saunders, W. C., Bouwes, N. & Jordan, C. Evaluating impacts using a BACI design, ratios, and a Bayesian approach with a focus on restoration. Environ. Monit. Assess. 188, 555 (2016).
PubMed Central  Article  Google Scholar 

77.
Murakami, D. spmoran (ver. 0.2.0): an R package for Moran eigenvector-based scalable spatial additive mixed modeling. Preprint at https://arxiv.org/abs/1703.04467v9 (2017).

78.
Murakami, D. & Griffith, D. A. Spatially varying coefficient modeling for large datasets: eliminating N from spatial regressions. Spat. Stat. 30, 39–64 (2019).
Article  Google Scholar 

79.
Murakami, D. & Griffith, D. A. Balancing spatial and non-spatial variation in varying coefficient modeling: a remedy for spurious correlation. Preprint at https://arxiv.org/abs/2005.09981 (2020).

80.
Walker, W. et al. Forest carbon in Amazonia: the unrecognized contribution of Indigenous territories and protected natural areas. Carbon Manag. 5, 479–485 (2014).
CAS  Article  Google Scholar 

81.
Robinson, E. J., Albers, H. J. & Busby, G. M. The impact of buffer zone size and management on illegal extraction, park protection, and enforcement. Ecol. Econ. 92, 96–103 (2013).
Article  Google Scholar 

82.
Koop, G. & Tole, L. Is there an environmental Kuznets curve for deforestation? J. Dev. Econ. 58, 231–244 (1999).
Article  Google Scholar 

83.
Barnes, M. D., Craigie, I. D., Dudley, N. & Hockings, M. Understanding local-scale drivers of biodiversity outcomes in terrestrial protected areas. Ann. NY Acad. Sci. 1399, 42–60 (2017).
PubMed  Article  PubMed Central  Google Scholar 

84.
Chamberlin, T. C. The method of multiple working hypotheses. Science 15, 92–96 (1890).
Google Scholar  More

Del Mar methane seep carbonate community analyses
The carbonate rock sample used for microbial community analyses was collected from Del Mar East Seep, Dive SO 177, (32°54.25456764 N, 117°46.9408327 W) at a depth of 1032.06 m. The site is in the northern portion of the San Diego Trough, about 50 km west of San Diego, California. Visible features of the seep include carbonate boulders and pavements colonized by orange and white bacterial mats, possible subsurface methane hydrate (large pits and craters), clam beds, and curtains of methane bubbles [20]. For a more in-depth discussion of the Del Mar Methane Seep, please see Grupe et al., (2015). A small chunk of the rock was sealed in a Mylar bag and shipped on dry ice to the University of Minnesota (Twin Cities). Upon arrival, the carbonate sample was temporarily stored at −80 °C until sampling for DNA extractions. All tools used for sampling biomass were autoclave-sterilized prior to use, and all work was performed in a LabConco A2 biological safety cabinet. Sterile aluminum foil was placed over ice packs to provide a cold and clean working environment. Three tubes were designated for the top of the rock, and three tubes were designated for the bottom of the rock. Biomass was scraped from the respective positions and placed in tubes. Top sample tubes had a biomass weight between 76.5–97.6 mg, and bottom sample tubes had a biomass weight between 30.7–59.3 mg. DNA was extracted from the seep samples using the ZymoBIOMICS DNA miniprep kit (Zymo, Irvine, CA). Nuclease-free water from the kit was processed alongside the samples as a negative control for iTag sequence analyses. The amplification of DNA and the generation iTag libraries of the V4 hypervariable region was performed by the University of Minnesota Genomics Center as previously described [21]. The samples and negative control were sequenced on ¼ of a lane of MiSeq for paired end 2×300 bp reads. Primers and adapters were removed with Cutadapt v. 2.10 [22]. The paired end reads were processed and assembled using DADA2 v1.16.0 [23]. The maximum expected error rate was 2 and reads detected as phiX were removed prior to error detection, merging of pairs and chimera detection. Taxonomic assignment was performed using the Silva database v. 138 [24]. Bioinformatic and statistical analyses was performed using tools in PhyloSeq 1.32.0 [25]. R package Decontam 1.8.0 [26] identified contaminating amplicon sequence variants (ASVs). Thus, ASVs of the genera Escherichia/Shigella, Haemophilus, Streptococcus, and the clade Chloroflexi S085 were bioinformatically removed from analyses. ASVs with statistically different abundances between the top and bottom surfaces were detected with DESeq2 v. 1.28.0.
Scanning electron microscopy of seep carbonates
The carbonate rock sample used for SEM imaging was collected from the Lasuen Knoll Seep (Dive SO 170, 33°23.57489996 N, 118°0.39814252 W) at a depth of 279.705 meters. A small chunk of carbonate rock with filtered bottom water and 25% vol/vol glutaraldehyde was stored at 4 °C, shipped on ice to the University of Minnesota (Twin Cities), and then stored at 4 °C. Several small pieces of the carbonate rock were broken off and placed in an 8-well plate. The carbonate pieces in the 8-well plate were rinsed, and subjected to an ethanol dehydration series as follows: 50% for 2 h, 70% over night, 80% for 15 min, 95% (x2) for 15 min, and 100% (x2) for 15 min. Carbonate samples were then subjected to critical point drying on a Tousimis Model 780 A Critical-Point-Dryer following standard procedures, followed by sputter coating 1–2 nm Iridium in a Leica ACE600 Sputter Coater. Lastly, carbonate samples were visualized at 1.0 kV on a Hitachi SU8230 Field Emission Gun Scanning Electron Microscope.
Continuous flow bioreactor experiments
To investigate the potential for sulfur-oxidizing bacteria to dissolve carbonate minerals, we used flow-through biofilm reactors (herein referred to as bioreactors: CDC; CBR 90–3 CDC Biofilm Reactor). Bioreactors contained eight polypropylene coupon holder rods that each accommodate three 12.7 mm diameter coupons. The lid and coupon holder rods were mounted in a 1 L glass vessel with side-arm discharge port at ~400 mL. A liquid medium (described below and Supplementary Table 1) was circulated through the bioreactor, while mixing was generated by a magnetic stir bar at 80 RPM. Sampling of the coupons was conducted aseptically by removing individual coupon holders and harvesting the coupons, while replacing the removed coupon holder with a sterile rubber plug. Experiments were 21 days in duration. On the final day of bioreactor experiments, one coupon was transferred to a treatment imaging flow cell (herein referred to as flow cell: BioSurface Technologies; Model FC 310 Treatment Imaging Flow Cell) in order to visualize and measure the in vivo pH of the biofilm. The flow cell is an autoclavable polycarbonate plastic cell with a recessed inner circle for installing coupons with biofilms grown in CDC Biofilm Reactors, and is equipped with barbed influent and effluent connectors, and topped with a coverslip for in vivo imaging.
The bioreactor medium consisted of two solutions, salt solution 1 & 2, and six additional components: vitamin solution, trace element solution, yeast extract, sodium bicarbonate, sodium thiosulfate, and sodium metasilicate. One 10 L carboy was prepared by autoclaving 5 L of salt solution 1 and one 10 L carboy was prepared by autoclaving 4.1 L (with an additional 839 mL ddH2O) of salt solution 2 at 121 °C for 80 min. The salt solutions were then aseptically pumped together, followed by the addition of sterile stocks of 10 mL 1000X vitamin solution, 10 mL 1000X trace element solution, 10 mL (1000 g/L) yeast extract, 11 mL 1 M NaHCO3, 10 mL 1 M Na2S2O3, and 20 mL 100 mM NaSiO3 through a 0.2 μM filter.
The final medium had a pH of ~7.85 and contained approximately 1.1 mM HCO3-, 1 mM S2O32-, and 9.14 mM Ca2+, with a saturation state with respect to aragonite (Ωaragonite) of ~0.5, similar to that measured near the sediment/water interface in certain seep environments where carbonate dissolution was observed to occur [7, 27]. The saturation state of aragonite was calculated using CO2SYS_v2.1 [28]. Ingredients for both salt solutions and the additional components added to the medium are presented in Supplementary Table 1. Measured saturation indices throughout bioreactor experiments are presented in Supplementary Table 2.
Mineral preparation for bioreactor experiments
Bulk mineral specimens of aragonite were obtained from D.J. Minerals (Butte, Montana). X-ray powder diffraction (XRD) of the aragonite sample, and peak matching indicates the mineral is nearly 100% aragonite. Samples used for the bioreactor were cored with a diamond coring bit on a drill press, cut to ~1–3 mm thickness and smoothed with a diamond saw. Coupons were then scored with a diamond tipped scribe pen to assign a reference code. After scoring, coupons were placed in a muffle furnace (NEY 2-525 Series II) at 500 °C for 4 h to remove residual organic matter, followed by weighing to 0.01 mg (CAHN 29 Automatic Electrobalance) and stored in a desiccator until being mounted in coupon holder rods.
Assembling the bioreactor
Weighed coupons were mounted in coupon holder rods, and fastened in place with plastic screws to avoid mineral chipping, as occurred with the stock metal screws. The bioreactor lid was equipped with a bubble trap on the media inflow port, a 0.45 μm airport, and a lure-lock covered in foil for inoculating via syringe. The entire bioreactor was autoclaved at 121 °C for 45 min.
Microbial inoculant
Celeribacter baekdonensis strain LH4 was used as the pure-culture inoculant for the bioreactor experiments. C. baekdonensis strain LH4 is a colorless, chemolithoheterotrophic, sulfur-oxidizing bacterium, belonging to the Rhodobacterales within the Alphaproteobacteria. It was isolated from marine sediments collected from a methane seep/brine pool at Green Canyon Block 233 [29] (Lat/Long: 27° 43.4392’ N, 91° 16.7638’ W, Depth: 648 m) in the Gulf of Mexico. Strain LH4 produces acid from thiosulfate oxidation, likely via the sox pathway (soxABCDXYZ) [30]. Strain LH4 also contains other genes related to sulfur oxidation, including four copies of sulfide quinone oxidoreductase ORFs and three copies that encode flavocytochrome c oxidases. Strain LH4 was grown up in the final medium, centrifuged at 10,000 G for 10 min, and washed 2X. OD590 for experiments was on average ~0.07 (~2.3•108 cells/mL based on cell counts and growth curves yielding the equation y = 3•109(x) + 2•107 where y equals cells/mL, and x equals OD590) with an inoculum volume of 20 mL.
Bioreactor experiments
Experiments were run to measure the dissolution of aragonite in a medium with a saturation state of 0.5 with respect to aragonite, similar to that measured near the sediment/water interface in seep environments [27]. Unlike in the sediments where AOM-generated alkalinity is high and carbonate actively precipitates, at the sediment/water interface previous studies have measured a notable drop in alkalinity [11, 31,32,33,34,35,36] and carbonate dissolution was observed to occur [7,8,9, 11, 27]. In total, six continuous -flow bioreactor experiments were run: three biotic experiments under identical conditions at 10 °C, pH ~7.85, 1.1 mM HCO3-, 1 mM S2O32-, and yeast extract; two uninoculated controls run under the same conditions as above to obtain an abiotic dissolution rate in undersaturated conditions; and one biotic experiment was run under identical conditions, minus the addition of S2O32- to elucidate the impact of heterotrophy on mineral stability.
10 L carboys containing the final medium and sterile bioreactors were assembled aseptically in a class A2 biological safety bench. Approximately 325 mL of the media was pumped into the bioreactor, followed by inoculating with 20 mL of pure culture and ~600 μL 1 M HCO3-. The luer-lock was removed and a 0.2 μM filter was added in its place. Bioreactors were then placed in a 10 °C refrigerator on a stir plate at 80 RPM.
Inoculated bioreactors were run in batch phase for 4 h to promote ample attachment of cells to aragonite coupon surfaces. Following the batch phase, the bioreactor was moved back to the clean hood where one coupon holder rod was removed and replaced with a sterile rubber stopper. Two coupons were placed in the incinerator at 500 °C for 4 h, followed by weighing for mass loss. The additional coupon was fixed in 4% PFA for 2 h, washed, and stained with DAPI [4,6-diamidino-2-phenylindole] for 30 min. Stained coupons were washed and mounted on coverslips (0.17 mm thickness) with DPX mountant for cell counting. After batch phase harvesting, the pump was turned on to the max flow rate (31 mL/min) and the first 100 mL out of the reactor was collected to measure pH, alkalinity, [Ca2+], and OD590. pH and alkalinity were measured using a Hanna Instruments total alkalinity mini titrator (HI-84431) and pH meter (HI 1131B) following the manufacturers protocol (Hanna Instruments, Woonsocket, RI, USA). [Ca2+] was measured using a Hach hardness test kit (product #2063900) following to the manufacturer’s protocol (Hach, Loveland, CO, USA). OD590 was measured on a Thermo Scientific Spectronic 20D + (Thermo Fisher Scientific, Waltham, MA, USA), where 1 mL of sterile medium was added to a cuvette to blank the spectrophotometer, and then 1 mL of the collected outflow was placed in a cuvette and measured. The pump was kept at the max flow rate for approximately 3–4 h after batch phase to remove planktonic cells and the additional HCO3- added during batch phase. Directly after flushing the bioreactor, water chemistry measurements were taken as they were above, and then the flow rate was dropped to 10 mL/min. Henceforth, the flow rate was only adjusted to keep the bulk fluid at approximately the pH and alkalinity of the starting conditions (max flow rate of 21 mL/min during week 3). To ensure consistency of the bulk fluid chemistry, pH and alkalinity measurements were taken twice daily (8–12 h intervals), just before changing carboys. Carboys were changed about every 8–12 h, depending on the current flow rate. To change carboys, the following steps were taken: (1) The pump was briefly turned off, and sterile aluminum foil was wrapped around the outflow tubing. (2) The bioreactor and currently connected carboy were moved to the biological safety bench on a roll cart. (3) Tubing that connected the carboy to the bioreactor was disconnected, and the male-port of the tubing still connected to the bioreactor was submerged in sterile 70% EtOH. (4) A new 10 L carboy containing the final medium was removed from the 10 °C refrigerator. (5) Aluminum foil covering the female-port of the tubing on the carboy was removed, and sprayed with sterile 70% EtOH. (6) Tubing from the bioreactor and carboy were then connected and placed on the roll cart, and moved back to the 10 °C refrigerator. (7) The pump was turned back on and flow resumed. One coupon holder rod was harvested every four days for 21 days, where two coupons were incinerated and weighed for mass loss and one coupon was fixed and stained for cell counting as described above. Surface area of aragonite coupons were estimated using the diameter of the coupon, the density of aragonite (2.93 g/cm3), and the mass of the coupon before and after bioreactor experiments. Density and mass were used to estimate the volume of the coupon, followed by deriving an estimated surface area. Using the mass of coupons before and after experiments, we calculated the total number of moles lost per unit area. Averaging these data out over time we were able to determine a dissolution rate in units of μmol CaCO3 • cm−2 • hr−1, which we then converted to μmol CaCO3 • cm−2 • yr−1 (8,760 h/year).
Confocal microscopy
We utilized the flow cell in conjunction with the ratiometric dye C-SNAFL-1 [5-[6]-carboxyseminaphthofluorescein] and Hoechst 33342 for measuring the in vivo pH of 21-day-old biofilms, and visualizing the biofilm, respectively. Coupons transferred to the flow cell were connected to a peristaltic pump at 3 mL/min until imaging occurred.
Spectrofluorometric assays
To calibrate C-SNAFL-1 photon excitations and to evaluate the dye’s potential use in measuring pH in biofilms, 1 mL aliquots of the final media with 10 mM HEPES buffer were adjusted to pH 5.0, 5.2, 5.6, 6.0, 6.4, 6.8, 7.2, 7.6, and 8.0. The final probe concentration was ~1.09 μM of C-SNAFL-1. Additional calibrations were performed with the addition of 4% (vol/vol) C. baekdonensis strain LH4.

$${mathrm{Ratio}} = left( {{mathrm{Ex}}_{488} – {mathrm{Ex}}_{488{mathrm{,bkgd}}}} right)/left( {{mathrm{Ex}}_{561} – {mathrm{Ex}}_{561{mathrm{,bkgd}}}} right)$$
(IV)

Confocal scanning laser microscopy
All biofilm images, pH and ratio calibrations in the flow cell were collected using a Nikon TiE inverted microscope equipped with an A1Rsi confocal scan head (Nikon). Hoechst 33342 fluorescence was excited with a 405 nm laser and collected between 425 and 475 nm using a regular PMT to visualize the biofilm. SNAFL emission was collected between 570 and 620 nm with a GaAsP PMT alternating excitation with the 488 and 561 nm lasers to obtain a ratiometric image. This image was used to calculate the pH across the field (below). 10X/0.45 PlanApo and 20X/0.75 PlanApo VC objectives were used to collect 512- by 512-bit resolution z-stacks. Depths in this paper are the distance from the substratum (base of biofilm) to the focal plane (distance furthest from the biofilm).
To visualize the biofilms, flow was suspended and coupons were stained with Hoechst 33342 (2 drops/mL medium) for 15 min and then flushed. C-SNAFL-1 was added to the medium at a final concentration of 20 μM. Flow resumed with C-SNAFL-1 in the medium and then was suspended for no more than 15 min to minimize the accumulation of acidic metabolites, which could artificially modify the local pH. Image planes were set up to visualize the coupon-attached biofilm and the bulk media adjacent to the coupon in the same image. Biofilms were excited at 405 nm, and emission was detected between 425 and 475 nm. pH microenvironment images were captured using dual excitation at 488 nm and 561 nm, and single emission at 600 nm. To determine the z-stack image parameters, we first excited the sample at 405 nm to visualize the base of the biofilm, and then switched to the dual excitation channels to determine the end of the focal plane (e.g. the coverslip; point furthest away from the biofilm in the z-dimension where the signal diminishes). After determining the z-distance to image, we imaged the biofilm at 405 nm, and then the pH of the microenvironment at 488/561 nm. The pH of the microenvironment was determined by calculating the ratio of excitation intensities (i.e., pixel values) between the two channels (Ex488 nm/ Ex561 nm). All intensities were determined using NIS Elements software.
Calibration standards for CSLM were performed in the flow cell with a sterile aragonite coupon and media with 20 μM C-SNAFL-1 and 10 mM HEPES buffer at pH 8.0, 7.8, 7.7, 7.55, 7.2, 7.0, and 6.8. Calibrations were performed with identical microscope settings as in the experiment above. Titration curves of pH versus ratios were calculated using NIS Elements software according to Eq. IV and were used to convert C-SNAFL-1 excitation ratios to pH values.
Image analysis
pH gradients within the biofilm and surrounding microenvironment were evaluated with NIS Elements software by using kymographs on the ratio images (Ex488 nm/ Ex561 nm). Each kymograph was ~500 μm in length (xy) and ~100 μm in depth (z). The resulting kymograph image had 4–5 regions of interest selected, including 2–3 spanning the biofilm into the bulk fluid directly above, and 2-3 in the bulk fluid adjacent to the biofilm. Regions of interest varied in z-depth depending on the thickness of the biofilm, though all were 13 μm in width (xy).
Estimation of CaCO3 dissolved from various seep locations
Seep locations
PDC1 – Dive, SO 162; Site, Point Dume; Feature, Chimney complex 1; Latitude, 33 56.45753693; Longitude, 118 50.70879299; Depth (m), 729.522
PDC2 – Dive, SO 164; Site, Point Dume; Feature, Chimney complex 2; Latitude, 33 56.50074983; Longitude, 118 50.64165898; Depth (m), 725.723
PDC3 – Dive, SO 164; Site, Point Dume; Feature, Chimney complex 3; Latitude, 33 56.45417243; Longitude, 118 50.70728425; Depth (m), 728.966
Determining total surface area
Each image taken by submersible at Point Dume had two scale bars generated at the time of imaging, one in the background and one in the foreground, with each representing 10 cm. These two scale bars, in conjunction with a designed macro in Fiji/ImageJ were used to interpolate the total surface area within each image. The macro is based on the equation y-y1 = m(x-x1). The slope, m, was determined using the equation ((y2-y1)/(x2-x1)), where y2 = 10 cm/(length in pixels of the foreground scale bar), y1 = 10 cm/(length in pixels of the distant scale bar), x2 = (y-coordinate of the foreground scale bar), x1 = (y-coordinate of the distant scale bar). Each image was imported into Fiji/ImageJ, and scaled to 1333 × 750 pixels. All scales were removed to put all measurements in pixel values. Centroid and integrated density were added to set measurements. Each image was duplicated and the duplicated images were changed to 32-bit images. The macro was then run on each respective image. A polygon was then drawn around the visible area within each image, and added to the ROI manager. The ROIs were added to their respective duplicate and measured. These data gave us the total visible surface area present in each image.
Determining surface area of exposed carbonates and carbonated-attached bacterial mats
Polygons were drawn around all exposed carbonates within the visible region. These ROIs were then added to their respective duplicated image, and the total area of exposed carbonates was measured. Carbonate ROIs that contained small regions in which bacterial mats were not visible had “negative” ROIs drawn around them. When applicable, these negative ROIs were subtracted from the total exposed carbonate measurement, resulting in the total mat-on-rock surface area measurement. It is likely that bacteria are also colonizing the surfaces on which no obvious mat can be observed, but removal of these areas from consideration was done to provide a conservative lower estimate of bacterial coverage.
Moles of carbonate rock dissolved was calculated as follows
Experimentally-derived carbonate dissolution rate = 1773.97 μmol CaCO3 • cm−2 • yr−1.
Apply dissolution rate to seep surface area covered by rock (e.g. Hydrate Ridge) to obtain the annual dissolution; 17.7397 mol CaCO3 • m−2 • yr−1 * 3.00E + 05 m2 = 5.32E + 06 mol CaCO3 • yr−1.
Apply 92.77% average carbonate-attached bacterial mat percent coverage (Supplementary Table 3) to annual dissolution; (5.32E + 06 mol CaCO3 • yr−1) • (0.9277) = 4.94E + 06 mol CaCO3 • yr−1.
We then applied these calculations to the remaining five seep sites in Table 1 to obtain the annual amount of carbonate dissolved per seep.
Table 1 Estimation of moles of carbon potentially released to the ocean/atmosphere system from the weathering of carbonate rocks via lithotrophic sulfur-oxidation at seep sites where carbonate coverage has been estimated.
Full size table More

The analytical model consists of (1) an adapted drag formulation for closely-packed cylinder arrays, including blockage and sheltering, and (2) a turbulent kinetic energy balance, necessary to quantify sheltering. The turbulence model builds on the formulation suggested by Nepf25 for vegetation canopies, and incorporates a turbulence production term by flow expansion, and an extended drag formulation in the wake production term. The steps to derive the equations are presented below.
Drag model
The drag forces experienced by an array of cylinders, per unit mass, can be calculated as:

$$begin{aligned} F_{d} = frac{1}{2}c_D a |U|U end{aligned}$$
(1)

where (c_D) is the drag coefficient of a single cylinder, which can be estimated using the empirical expression of White30, given by:

$$begin{aligned} c_D = 1 + 10Re^{-2/3} end{aligned}$$
(2)

where Re is the Reynolds number based on the cylinder diameter and the depth-averaged local flow velocity U. a is the projected plant area per unit volume, defined by Nepf25 as:

$$begin{aligned} a = frac{d h}{h s^2} = frac{d}{s^2} end{aligned}$$
(3)

with d being the cylinder diameter, s the spacing between cylinders, and h the water depth.
The main unknown in Eq. (1) is the local flow velocity U. If a cylinder array is sufficiently sparse, the local flow velocity could be assumed equal to the depth-averaged incoming flow velocity, (U_{infty }), either measured or calculated with a free surface flow model. For denser configurations, the velocity will change as the flow propagates through the array due to (1) flow acceleration between the elements (blockage), and (2) flow deceleration due to the sheltering effects of upstream rows of cylinders. Both effects are illustrated in Fig. 1c. The changes in flow velocity are included by multiplying (U_{infty }) by a blockage factor, (f_b), and a sheltering factor, (f_s):

$$begin{aligned} U = f_b f_s U_{infty } end{aligned}$$
(4)

Inserting both factors in the expression for the drag force results in Eq. (5):

$$begin{aligned} F_{d} = frac{1}{2}c_D a |U|U = frac{1}{2}c_D a f_b^2 f_s^2 |U_{infty }|U_{infty } = frac{1}{2} c_{D,b} a |U_{infty }|U_{infty } end{aligned}$$
(5)

where the changes in velocity have been incorporated in the bulk drag coefficient, (c_{D,b} = c_D f_b^2 f_s^2). This expression provides a direct relationship between the drag coefficient of a single cylinder, (c_D), and bulk drag coefficients (c_{D,b}) measured for cylinder arrays in laboratory experiments. Predicting the drag force thus depends on determining the values of (f_b) and (f_s).
The blockage factor (f_b) can be estimated based on mass conservation through a row of cylinders11, considering that the velocity will increase as the same flow discharge travels through the smaller section between the elements:

$$begin{aligned} U_{infty } A = U_c A_c = f_b U_{infty } A_c end{aligned}$$
(6)

where the total frontal area is (A = h s_y), and (s_y) is the distance between cylinders perpendicular to the flow, center-to-center (see Fig. 1). Subtracting the frontal area of the cylinders from the total area gives the available flow area, (A_c):

$$begin{aligned} A_c = h s_y – h D = h (s_y-d) end{aligned}$$
(7)

Here we are assuming that the water depth is the same just upstream and in between the cylinders. Solving for (f_b) in Eq. (6) results in Eq. (8), see also Etminan et al.11:

$$begin{aligned} f_b = frac{h s_y}{ h (s_y-d)} = frac{1}{1-d/s_y} end{aligned}$$
(8)

The sheltering factor (f_s) can be estimated from the wake flow model by Eames et al.26, which predicts the velocity deficit behind a cylinder as a function of the distance downstream of the cylinder, (s_x), the cylinder diameter, the local turbulent intensity (I_t), and the drag coefficient:

$$begin{aligned} frac{U_{infty }-U_{w}}{U_{infty }} = frac{c_D d}{2sqrt{2 pi } I_t s_x} end{aligned}$$
(9)

where (U_{w}) is the velocity in the cylinder wake, (U_{infty }) is the incoming flow velocity, and (I_t) is the meant turbulent intensity, defined as (I_t = sqrt{k}/U_{infty })21,25. k represents the turbulent kinetic energy per unit mass, with (k = 1/2(overline{u’^{2}} + overline{v’^{2}} + overline{w’^{2}})), where (u’), (v’), and (w’) are the instantaneous velocity fluctuations in the streamwise, lateral, and vertical direction respectively, and where the overbar denotes time averaging. The turbulent velocity fluctuations are defined as the difference between the instantaneous velocities and their mean value over a measurement period. Here we consider the depth-averaged value of the turbulent intensity, in view of the uniformity of the turbulent properties over the vertical observed inside emergent arrays25.
Equation (9) was developed assuming turbulent flow. Viscous effects decrease the velocity deficit26, with the reduction factor being given by:

$$begin{aligned} f_{Re} = sqrt{frac{Re}{Re_{t}}} end{aligned}$$
(10)

where (Re_{t}) is the lowest Reynolds number corresponding to fully turbulent wake flow. Laminar effects are included in the wake flow model by multiplying the velocity deficit of Eq. (9) by the reduction factor (f_{Re}) for (Re < Re_t), where the the turbulent Reynolds number is assumed equal to (Re_t = 1,000). This value is based on the observation that although a wake starts becoming turbulent at (Re_{t} sim 200), drag coefficient measurements usually become constant at Reynolds numbers beyond (Re_{t} sim 1000), e.g. as shown in Figure 2.7 of Sumer and Fredsoe13. The influence of varying (Re_{t}) on the model results is investigated in “Results and discussion” section. Defining the sheltering factor as (f_s = frac{U_{w}}{U_{infty }}), and including (f_{Re}) and the bulk drag coefficient in the definition of the velocity deficit results in Eq. (11): $$begin{aligned} f_s = frac{U_{w}}{U_{infty }} = 1-f_{Re}frac{c_{D,b} d}{2sqrt{2 pi } I_t s_x} = 1-f_{Re}frac{c_{D,b} d}{2sqrt{2 pi } (sqrt{k}/U_{infty }) s_x} end{aligned}$$ (11) Equation (9) also assumes that the downstream cylinder is placed inside the ballistic spreading region of the wake. The ballistic regime occurs for a downstream distance (s_x < L/It), where L is the integral length-scale of turbulence, and it is characterized by a rapidly decaying velocity deficit, and by a linear increase of the wake width with downstream distance. Inside the cylinder arrays, the length scale development is limited by the downstream spacing, resulting in (L < s_x). Considering that turbulent intensity measurements of Jansen29 varied between (I_t) = 0 and 0.8 inside cylinder arrays with n = 0.64–0.9, this would result in (L < s_x/It). This is a reasonable general assumption for the bamboo structures, since their porosity varies in a similar range. If the poles were sparsely placed, there would be a transition from ballistic to diffusive spreading of the wake. Eames et al.26 also developed expressions for turbulent flow under the diffusive regime, which could be used in place of Eq. (9). In the opposite case of very high pole densities, there may be a point where the elements are so closely-packed that vortex shedding is inhibited by the presence of the neighboring cylinders. Considering an analogy with a cylinder placed close to a solid boundary, vortex shedding would not take place for spanwise spacings smaller than (s_y/d < 1.3)13, causing a decrease of the drag coefficient that would not be reproduced by the expression of White30. The application of the present model is thus restricted to (s_y/d > 1.3).
Balance of turbulent kinetic energy
Application of Eq. (11) requires predicting the turbulent kinetic energy. This is calculated by expanding the model developed by Nepf25, based on a balance between turbulence production and dissipation:

$$begin{aligned} P_w sim epsilon end{aligned}$$
(12)

where (P_w) is the turbulent production rate and (epsilon) is the dissipation rate. For a dense cylinder array, k is produced by (1) generation in the wakes of the cylinders25, and (2) shear production by the jets formed between the elements28. The total turbulence production term, (P_w), consequently has two parts:

$$begin{aligned} P_w = P_{w1}+P_{w2} end{aligned}$$
(13)

We assume that for dense cylinder arrays these two terms are much higher than turbulence production by shear at the bed, based on observations by Nepf25 for sparse arrays. This assumption is further tested in “Results and discussion” section.
The first term in Eq. (13), (P_{w1}), represents turbulence production at the wakes, and can be estimated as the work done by the drag force times the local flow velocity:

$$begin{aligned} P_{w1} = frac{1}{2}c_D a |U|U^2 = frac{1}{2}c_D a f_b^3 f_s^3 |U_{infty }|U_{infty }^2 end{aligned}$$
(14)

The second term, (P_{w2}), represents turbulence generation due to flow expansion28, and can be estimated from the Reynolds shear stresses:

$$begin{aligned} P_{w2} = overline{ u’ v’} frac{partial u }{partial y} end{aligned}$$
(15)

where the overbar denotes time averaging. The loss in mean kinetic energy (E_c) due to flow expansion is equal to:

$$begin{aligned} Delta E_c = frac{1}{2} U_{infty }^2 left( left( frac{A}{A_c}right) ^{2}-1 right) = frac{1}{2} left( f_b^{2}-1 right) U_{infty }^2 end{aligned}$$
(16)

where the energy loss due to flow expansion, (Delta E_c), is modelled using the Carnot losses. Assuming that the mean kinetic energy is transformed into turbulent kinetic energy (E_t), and assuming isotropic turbulence, gives Eq. (17):

$$begin{aligned} frac{1}{2} left( f_b^{2}-1 right) U_{infty }^2 = frac{3}{2}overline{ u’ u’} end{aligned}$$
(17)

Equation (17) enables expressing the normal Reynolds stress as a function of the incoming flow velocities and the blockage factor:

$$begin{aligned} overline{ u’ u’} = frac{1}{3} left( f_b^{2}-1 right) U_{infty }^2 end{aligned}$$
(18)

The Reynolds shear stress is estimated as (overline{ u’ v’} = Roverline{ u’ u’}), where the correlation factor R was given a constant value of 0.4 based on observations of Nezu and Nakagawa31. This value was derived for open channel flow conditions and is assumed acceptable as a first approximation, but it could vary inside a cylinder array. This is explored further in “Results and discussion” section.
The velocity gradient is estimated from the velocity difference between the side of the cylinders (dominated by blockage) and the wake of the cylinders (dominated by sheltering) resulting in Eq. (19):

$$begin{aligned} frac{partial u }{partial y} approx frac{U_{infty }(f_b-f_s)}{frac{1}{2} s_y} end{aligned}$$
(19)

Substitution into Eq. (15) gives Eq. (20):

$$begin{aligned} P_{w2} = frac{2}{3} R (f_b-f_s)(f_b^{2}-1)frac{U_{infty }^3}{s_y} end{aligned}$$
(20)

The dissipation term, (epsilon), is estimated as:

$$begin{aligned} epsilon sim k^{3/2} l^{-1} end{aligned}$$
(21)

The characteristic turbulent length scale l is limited by the surface-to-surface separation between the elements in the flow direction, (l = min(|s_x-d|, d)). This differs from the expression developed by Nepf25, who used the diameter as representative of the size of the eddies. We assume that in closely-packed cylinder arrays the spacing between cylinders may be smaller than the diameter, (|s_x-d| < d), which would limit turbulence development. The maximum value of l is set equal to the cylinder diameter. Here we also assume that for the dense cylinder arrangements, the spacing between cylinders is considerably smaller than the water depth, hence turbulence generated by bed friction is negligible. Balancing the production and dissipation of turbulent kinetic energy results in Eq. (22): $$begin{aligned} frac{k^{3/2}}{l} sim |U_{infty }|U_{infty }^2left( c_D a f_b^3 f_s^3 + frac{ 4R}{3s_y}(f_b^{2}-1)(f_b-f_s)right) end{aligned}$$ (22) Taking the cubic root at both sides and introducing the scale factor (alpha _1) gives Eq. (23): $$begin{aligned} frac{sqrt{k}}{U_{infty }} = alpha _1left( c_D f_b^3 f_s^3 a l + frac{4}{3}R(f_b^{2}-1)(f_b-f_s)frac{ l}{s_y}right) ^{1/3} end{aligned}$$ (23) Where (alpha _1) is a coefficient of ({mathcal {O}}(1)), which is given a default value of (alpha _1 = 1). The sensitivity of the model to different (alpha _1) and R values is explored in “Results and discussion” section. k can be calculated by solving Eq. (23) iteratively, using the incoming upstream velocity (U_{infty }) and the geometric characteristics of the structure, (s_y, s_x, d) and a, as an input. This enables determining the sheltering factor, (f_s = U_{w}/U_{infty }) from Eq. (11). The blockage factor (f_b=(1-d/s_y)^{-1}) can also be calculated from the geometric properties of each configuration. Both coefficients can be then combined to predict the bulk drag coefficient, with (c_D,_{b} =c_D(f_s)^2(f_b)^2). Deriving (c_D,_{b}) with the present approach relies on the assumption that the changes in water depth through the structure are small. This is a reasonable assumption given the short length of the bamboo structures in the streamwise direction, which varies between 0.7 and 1.5 m (see Fig. 1b). Longer structures that experience non-negligible changes in flow depth and velocity should be discretized, and the bulk drag coefficient should be calculated separately for the different sections. The model assumptions are discussed further in the following section. More

Seeds collection and germination
We collected T. sebifera seeds by hand from populations in both the introduced (US—16 populations in total) and native (China—14 populations in total) ranges (for details see Table S1). At each population, we haphazardly selected 5–10 trees, and harvested thousands of seeds from each tree. In the laboratory, we removed the waxy coats around these seeds by hand after immersing them in a mixture of water and laundry detergent (10 g/L) for 24 h [29]. Then, we rinsed, surface sterilized (10% bleach), and dried them. In order to improve germination, we put these seeds in wet sand and stored them in the refrigerator (4 °C) for at least 30 days. In spring, we sowed these seeds in greenhouse trays (50 holes/tray) which were filled with sterilized (autoclave at 121 °C for 30 min) commercial potting soil, and then kept them in an open-sided greenhouse at Henan University in Kaifeng, Henan, China (34°49′13′′ N, 114°18′18′′ E) or unheated greenhouse at Rice University, Houston, TX USA (29°43′08′′ N 95°24′11′′ W). After seeds germinated and seedlings reached the 4 true leaf stage, we selected similar size seedlings to conduct the following experiments.
Common garden experiment—differences in AM fungal colonization and plant growth
To investigate the differences in AM fungal colonization and growth between plants from introduced (US) and native populations (CH), we carried out a common garden experiment at Henan University. We collected soil in a corn field, which includes most common AM fungal species based on previous reports [33, 34]. It was a sandy soil with total nitrogen and total phosphorus of 1.9 g/kg (DW) and 0.6 g/kg (DW), respectively, and pH of 7.68. We removed surface litter before collecting topsoil (10–15 cm depth) and then combined equal parts of soil and fine sand in 132 pots (21 cm × 16 cm, ~3 kg of soil mix each) after they were passed through a 1-cm mesh screen. We planted seedlings from 22 populations (12 native and 10 introduced populations, 6 seedlings of each population, Table S1) individually in these prepared pots and placed them in the open-sided greenhouse. We protected them from herbivores with nylon mesh (16 openings/cm) cages during the experiment. After 60 and 90 days, we harvested 3 seedlings from each population as 3 reps each time and carefully washed their whole roots from the soil. We collected ~30 fresh fine roots ( >1 cm/segment) from each plant root to measure AM fungal colonization. In brief, we cleared (in 10% KOH), bleached, acidified, and stained (trypan blue) these samples then slide mounted 30 one cm long pieces of fine root for each plant [7]. AM fungal colonization was determined by the gridline intersect method with 300 intersection points per plant [35]. We dried and weighed the roots and shoots.
Collection of root exudates and flavonoids analysis for root exudates
Our previous study found higher concentrations of flavonoids but lower concentrations of tannins in roots of introduced populations of T. sebifera than in native populations [17] with quercetin and quercitrin being the main flavonoids [28, 30]. In our pilot experiment, we only detected quercetin and quercitrin in root exudates but no other flavonoids. Therefore, in this study we focused on quercetin and quercitrin in root exudates and their functions. We determined their amounts in root exudates from native (China) and introduced (US) populations at Henan University. We filled 132 glass beakers (1000 ml) with Hoagland’s solution [36] and covered the opening with a foam board with a hole in its center. We took 6 seedlings from each of 22 populations (12 native, 10 introduced, Table S1) and carefully washed the soil from their roots with tap water, then transplanted them individually into the beakers (1 seedling per beaker) and fixed them with a sponge. Because of mortality, only 80 plants of 17 populations (9 native, 8 introduced) survived until exudate collection. The odds of a plant dying did not depend on population origin (F1,20 = 3.7, P = 0.0679) or population (Z = 1.3, P = 0.0937). We checked these glass beakers and filled them with Hoagland’s solution every day.
After these plants grew for 57 or 87 days in an open-sided greenhouse with a typical temperature range of 18 °C (night) to 28 °C (day) and 13–14 h of natural daylight, we put DI water into these beakers instead of Hoagland’s solution to minimize the effects of environments on root exudates. Three days later (i.e., at 60 and 90 days) these plants were harvested to obtain their dry root mass. The root exudates were dried at 40 °C under vacuum by rotary evaporators. Then we extracted the chemicals from these concentrates in 3 ml of methanol solution with 0.4% phosphoric acid water (48:52, v:v) and filtered them through 0.22 μm hydrophobic membranes. The concentrations of quercetin and quercitrin were assessed by high-performance liquid chromatography [30]. In brief, 20 μl of extract was injected into an HPLC with a ZORBAX Eclipse C18 column (4.6 × 250 mm, 5 μm; Agilent, Santa Clara, CA, USA) with the following flow: 1.0 mL min−1 with a 100% methanol (B) and 0.4% phosphoric acid in water (A) as the mobile phase. The gradient was as follows: 0–10 min 52:48 (A:B); 10–24 min 48:52 (A:B). UV absorbance spectra were recorded at 254 nm. The concentrations of flavonoid compounds were calculated and standardized by peak areas of standards of known concentrations.
Root exudate addition experiment—effects of different populations on AM fungal colonization
In order to investigate the role that root exudates play in the interactions between AM fungi and plants, we conducted an experiment in which exudates were collected from plants in liquid (donor) and applied to the soils of other plants (target). The exudate donor plants were grown in 1080 (two venues: 540 seedlings at Rice University and 540 seedlings at Henan University) containers, each with 1000 ml of Hoagland’s solution, that each had a foam board top with a hole and a bottom drain tube that could be regulated. At each venue, we washed the soil from ~500 sets of plants (US = 465, China = 504) from native (8 populations for venue US and 7 populations for venue CH) or introduced (13 populations for venue US and 12 populations for venue CH) populations and secured them (3 plants per container) in the containers using sponges (details in Table S1). The remaining containers were left as plant-free controls. We started the application experiment after 7 days.
For exudate target plants, we collected the soil from different sites in the introduced or native ranges (See Table S1). At each site, we collected soil under the canopy of a T. sebifera tree (Home soil) and that more than 3 meters away from the canopy of a T. sebifera tree (Away soil). We collected the topsoil to a depth of 15 cm after removing the surface litter, air dried them, and screened them (1 cm mesh). These soils were mixed with vermiculite (1:2 volume). Then we used these mixes to fill 1080 pots (15 cm × 12 cm; 540/venue). Each pot in China received a mixed soil from a site in China and each pot in US received soil from a single small area within a site in the US. We transplanted a seedling from a native (12 populations for venue US and 3 populations for venue CH, See Table S1) or introduced (13 populations for venue US and 5 populations for venue CH, See the Table S1) population into each pot (270 of each per venue). We randomly assigned a target plant to each set of donor plants or water only controls.
Every 4 days we changed the Hoagland’s solution to DI water for 3 days to collect root exudates from donor plants. Then we applied this water solution from a donor set to its target plant. After 70 days, we harvested the target seedlings, kept a fine root sample for AM fungal colonization determination, then dried and weighed leaves, stems, and roots.
Chemical addition experiment—quercetin and quercitrin effects on AM fungal colonization
We transplanted 391 seedlings from 8 native populations (CH) and 9 introduced populations (US) into 391 pots with field soil (1.3 kg/pot) in nylon mesh cages at Henan University. To test the effect of quercetin and quercitrin on AM fungal colonization, we prepared solution of quercetin or quercitrin in acetone (10 mg/mL) (acetone did not affect AM fungal colonization based on our preliminary experiment). Then these solutions were diluted in water to 2 concentrations (1 mg/L and 10 mg/L) based on the result of chemical analyses of root exudates and the 0.1% of acetone in water as control. We watered 15 ml of solution (5 reps per population) or water (3 reps per population) around the base of seedling stems every 3 days (16 times in total). Four plants died (3 in quercitrin application treatment, 1 in quercetin application treatment). After 70 days, we collected seedlings by cutting at ground level and collected fine roots to test AM fungal colonization.
Activated carbon experiment—AM fungal colonization with inactivated chemicals
In order to verify the chemicals in root exudates play a key role in the relationship between AM fungi and plant roots, we conducted an experiment at Henan University with activated carbon (AC) addition to block bioactivity of root exudate chemicals. We filled plastic pots in mesh cages at Henan University with either 1.3 kg of field soil (N = 78) or field soils amended with activated carbon (N = 78, Sinopharm Chemical Reagent Co., Ltd, Beijing, China) added as 1:500 v:v. We transplanted seedlings from 13 populations (6 native and 7 introduced, Table S1) into the pots with 6 seedlings for each population. Eighteen seedlings died during this experiment (12 seedlings from AC treatment, 6 seedlings from control). After 70 days, we harvested plants and used a fine root sample to determine AM fungal colonization.
Field survey of AM fungal assemblages
We collected rhizosphere soil from 3 sites in China (Dawu, Hubei, 31°28′N, 114°16′E; Wuhan, Hubei, 30°32′N, 114°25′E; Guilin, Guangxi, 25°04′N, 110°18′E) for AM fungal species identification via high-throughput sequencing. At each of these sites, we selected 3 T. sebifera trees per site and dug the soil close to the tree trunk until its root branch was found. We collected soils from 3 roots per plant. We removed the bulk soil from these roots by shaking, and then collected the soil remaining on these roots using brushes (1 new brush per tree). The rhizosphere soils on the roots from same tree were mixed fully. About 5 g of fresh rhizosphere soil from one tree was collected and stored in dry ice and ultra-low temperature freezer (−80 °C) until they were used to test the AM fungi abundance based on high-throughput sequencing [37, 38].
For DNA extraction, microbial DNA was extracted from the prepared samples (0.25 g soil per sample) using the E.Z.N.A.® soil DNA Kit (Omega Bio-tek, Norcross, GA, U.S.) according to the manufacturer’s protocols. The DNA concentration and purification were determined by NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, USA), and DNA quality was checked by 1% agarose gel electrophoresis [39].
For the PCR amplification, nested PCR was conducted to amplify the AM fungi 18S rRNA. The primer pairs AML1 (5′-ATCAACTTTCGATGGTAGGATAGA-3′) and AML1 (5′-GAACCCAAACACTTTGGTTTCC-3′) were used in the first run. The primer pairs AMV4.5NF (5′-AAGCTCGTAGTTGAATTTCG-3′) and AMDGR (5′-CCCAACTATCCCTATTAATCAT-3′) were used in the second run in the thermocycler PCR system (GeneAmp 9700, ABI, USA). The PCR reactions were conducted using the program according to Xiao et al. [39].
For each sample, purified amplicons were pooled in equimolar and paired-end sequenced (2 × 300) on an Illumina MiSeq platform (Illumina, San Diego, USA) according to the standard protocols of Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The raw fastq files were quality-filtered by Trimmomatic and merged by FLASH with the following criteria: (i) the reads were truncated at any site receiving an average quality score More

Experimental design
The Winchmore Irrigation Research Station is in the centre of the Canterbury plains, the largest area of flat land in New Zealand (43.787° S, 171.795° E; Fig. 1). It is at an altitude of 160 m above sea level, a mean annual temperature of 12 °C, and has an annual rainfall of 745 mm (range 491–949 mm)20. The soil is a Lismore stony silt loam classified as an Orthic Brown soil in the New Zealand soil classification and as an Udic Ustochrept in USDA soil classification21. Flood irrigation, known locally as border-check/dyke irrigation, was installed at the site in 1947. However, the two long-term trials, hereafter known as the fertiliser and irrigation trials, were established in 1952 and 1949, respectively.
Fig. 1

Location of Winchmore within the Canterbury region (coloured green) and the layout of the long-term fertiliser and irrigation trials over time.

Full size image

Full details of the setup of the fertiliser and irrigation trials between 1949–1951, including the political rationale for the trial, its statistical design, cultivation dates, sowing rates of perennial ryegrass (Lolium spp) and white clover (Trifolium repens) and initial fertiliser and irrigation treatments are available elsewhere20.
The fertiliser trial has 20 border check irrigation bays divided into five treatments each with four replicates set out in a randomised block design (Fig. 1). From 1952/53 to 1957/1958 treatments were: nil (no P applied), 188, 376 (annually and split P applications), and 564 kg SPP ha−1. All P applications occurred annually in autumn except for the 376 kg SSP ha−1 treatment which had two treatments divided into an annual autumn application and split applications in between autumn and spring. From 1958/59 to 1979–80 the nil and 188 and 376 (split autumn and spring application) SSP treatments were unaltered, while P applications were stopped to the annual 376 and 564 SSP treatments. In 1972, 4.4 t/ha of lime (caclium carbonate) was applied to all treatments22. From 1980 onwards the nil, and 188 SSP treatments and the 376 SSP treatment, now receiving winter fertiliser applications, were joined by a treatment applying 250 SSP ha−1 in winter to the previous 376 SSP treatment and a Sechura rock phosphate treatment applying 22 kg P ha−1 in winter to the former 576 SSP treatment.
Each irrigation bay was fenced off, 0.09 ha in size and grazed by separate mobs of sheep at 6, 11, and 17 stock units (with one stock unit equivalent to one ewe at 54 kg live-weight) per replicate for the nil, 188 SSP, and 376 SSP treatments, respectively. This separation prevented carry-over of dung P and other nutrients and contaminants between treatments. No grazing occurred in winter. Flood irrigation was applied when soil moisture content (w w−1) fell below 15% (0–100 mm depth). This occurred on-average 4.3 times per year.
The irrigation trial had 24 irrigation bays (each 0.09 ha in size) which had lime applied to the whole trial in 1948 (5 t ha−1) and 1965 (1.9 t ha−1) to maintain soil pH at 5.5–6.0. From 1951 to 1954 treatments were SSP applied at 250 kg ha−1 in autumn annually and either four replicates of dryland, or five replicates of irrigation applied at one, two, three, six-weekly intervals or at three-weekly intervals in alternate seasons. From 1953/54 to 1956/57 the weekly and two-weekly treatments were replaced by irrigation when soil moisture in the top 100-mm of soil reach 50 and 0% available soil moisture (asm), respectively. In 1958 the irrigation trial was cultivated with a rotary hoe and grubber, 140 kg SSP ha−1 applied and the site re-sown in ryegrass and white clover. From 1958/59–2007 the site had the same blanket application of SSP and four replicates of dryland, while a completely randomised design was used to impose five replicates of four treatments (Fig. 1) that looked at irrigation applied when soil moisture in the top 100-mm of soil reach 10, 15 and 20% (equivalent to 50% asm with 0% asm being wilting point) and irrigation on a 21-day interval. The need for irrigation to the irrigation and fertiliser trials was informed by soil moisture measured weekly by technical staff using a mixture of gravimetric analyses (1950–1985), neutron probe (1985–1990) and time-domain reflectometer (1990-onwards). Irrigation was applied at a rate of 100 mm per event20.
Except for winter, when no grazing occurred, each treatment was rotationally grazed by a separate flock of sheep with 6 and 18 stock units per replicate for the dryland and 20% v/v treatments, respectively.
The irrigation trial finished in October 2007 although the P fertiliser regime continued. All irrigated treatments shifted to the same three weekly schedule as the long-term Fertiliser trial. The dryland treatment remained unirrigated. The Winchmore farm was converted into a commercial irrigated farm operation and sold in 2018. The fertiliser trial was also sold but with a covenant ensuring it continues to operate as per normal except that irrigation from 2018 onwards is now applied by spray irrigation with the aim of ensuring soil moisture is maintained above 90% of field capacity. Since January 2019 there are daily soil moisture meter records from a moisture meter installed into one of the control plots. Soil moisture, rainfall and irrigation are recorded.
The production of the Winchmore trials data records23 involved a three-step process (Fig. 2).
Fig. 2

Flowchart of the steps involved in sampling, analysis, collation and curation and data analysis and processing of the databases from the Winchmore Trials. Note that blue and orange boxes are sub tasks associated with each step and resulting outputs, respectively.

Full size image

Step 1: Soil and pasture sampling
Pasture production was measured from two exclusion cages (3.25 m long × 0.6 m wide) per plot24. Areas within each cage were trimmed to 25 mm above ground level and left for a standard grazing interval for that time of year. Following grazing a lawnmower was used to harvest a 0.40 m wide strip in the middle of each enclosure to 25 mm above ground level. The wet weight was determined, and a sub-sample taken to determine dry matter content with a separate sample manually dissected into grass, clover and weeds. All surplus mown herbage was returned to the plot. Approximately 9–10 cuts were made annually. A composite soil sample of 10 cores (2.5 cm diameter and 7.5 cm deep) was collected from each plot. These were collected four times annually (July, prior to fertiliser application, and October, January and April), using established best practices24,25. In 2009 soil samples were also collected from the 0–75, 75–150, 150–250, 250–500, 500–750, and 750–1000 mm depths on both trials17. During 2018, prior to cultivation, soil on the unirrigated, 10 and 20% soil moisture treatments of the irrigation trial were sampled at 0–150, 150–250, 250–500, 500–750, 750–1000, 1000–1500, and 1500–2000 mm depths. The top 250 mm of these samplings were collected by hand using an auger, but deeper depths were excavated via a mechanical digger. Representative sub-samples were taken from each depth. Annual samplings were crushed, dried and sieved More

1.
Gobler CJ, Doherty OM, Hattenrath-Lehmann TK, Griffith AW, Kang R, Litaker W. Ocean warming has expanded niche of toxic algae. Proc Natl Acad Sci USA. 2017;114:4975–80.
CAS  PubMed  Article  Google Scholar 
2.
Olivieri ET. Colonization, adaptations and temporal changes in diversity and biomass of a phytoplankton community in upwelled water off the Cape Peninsula, South Africa, in December 1979. South Afr J Mar Sci. 1983;1:77–109.
Article  Google Scholar 

3.
Irwin AJ, Zoe V, Finkel ZV, Müller-Karger FE, Troccoli, Ghinaglia L. Phytoplankton adapt to changing ocean environment. Proc Natl Acad Sci USA. 2015;112:5762–6.
CAS  PubMed  Article  Google Scholar 

4.
Chivers W, Walne A, Hays G. Mismatch between marine plankton range movements and the velocity of climate change. Nat Commun. 2017;8:14434.
CAS  PubMed  PubMed Central  Article  Google Scholar 

5.
Gisselson L, Granéli E, Pallon J. Variation in cellular nutrient status within a population of Dinophysis norvegica (Dinophyceae) growing in situ: single – cell elemental analysis by use of a nuclear microprobe. Limnol Oceanogr. 2001;5. https://doi.org/10.4319/lo.2001.46.5.1237.

6.
Ackermann M. A functional perspective on phenotypic heterogeneity in microorganisms. Nat Rev Microbiol. 2015;13:497–508. https://doi.org/10.1038/nrmicro3491.
CAS  Article  PubMed  Google Scholar 

7.
Núñez-Milland DR, Baines SB, Vogt S, Twining BS. Quantification of phosphorus in single cells using synchrotron X-ray fluorescence. J Synchrotron Radiat. 2010;17:560–6.
PubMed  PubMed Central  Article  CAS  Google Scholar 

8.
Berthelot H, Duhamel S, L’Helguen S, Maguer JF, Wang S, Cetinić I, et al. NanoSIMS single cell analyses reveal the contrasting nitrogen sources for small phytoplankton. ISME J. 2019;13:651–62. https://doi.org/10.1038/s41396-018-0285-8.
CAS  Article  PubMed  Google Scholar 

9.
Štrojsová A, Vrba J. Short-term variation in extracellular phosphatase activity: possible limitations for diagnosis of nutrient status in particular algal populations. Aquat Ecol. 2009;43:19–25.
Article  CAS  Google Scholar 

10.
O’Donnell DR, Hamman CR, Johnson EC, Kremer CT, Klausmeier CA, Litchman E. Rapid thermal adaptation in a marine diatom reveals constraints and trade-offs. Glob Change Biol. 2018;24:4554–65.
Article  Google Scholar 

11.
Jin P, Agustí S. Fast adaptation of tropical diatoms to increased warming with trade-offs. Sci Rep. 2018;8:17771. https://doi.org/10.1038/s41598-018-36091-y.
CAS  Article  PubMed  PubMed Central  Google Scholar 

12.
Thomas CD, Cameron A, Green RE, Bakkenes M, Beaumont LJ, Collingham YC, et al. Extinction risk from climate change. Nature. 2004;427:145–8.
CAS  PubMed  Article  Google Scholar 

13.
Urban MC. Accelerating extinction risk from climate change. Science. 2015;348:571–3.
CAS  PubMed  Article  Google Scholar 

14.
Kottuparambil S, Jin P, Agusti S. Adaptation of Red Sea Phytoplankton to experimental warming increases their tolerance to toxic metal exposure. Front Environ Sci. 2019;7. https://doi.org/10.3389/fenvs.2019.00125.

15.
Flores-Moya A, Costas E, Lopez-Rodas V. Roles of adaptation, chance and history in the evolution of the dinoflagellate Prorocentrum triestinum. Naturwissenschaften. 2008;95:697–703.
CAS  PubMed  Article  Google Scholar 

16.
Flores-Moya A, Rouco M, García-Sánchez MJ, García-Balboa C, González R, Costas E, et al. Effects of adaptation, chance, and history on the evolution of the toxic dinoflagellate Alexandrium minutum under selection of increased temperature and acidification. Ecol Evol. 2012;2:1251–9. https://doi.org/10.1002/ece3.198.
Article  PubMed  PubMed Central  Google Scholar 

17.
Martiny AC, Ustick LA, Garcia C, Lomas MW. Genomic adaptation of marine phytoplankton populations regulates phosphate uptake. Limnol Oceanogr. 2019. https://doi.org/10.1002/lno.11252.

18.
Ribeiro S, Berge T, Lundholm N, Andersen TJ, Abrantes F, Ellegaard M. Phytoplankton growth after a century of dormancy illuminates past resilience to catastrophic darkness. Nat Commun. 2011;2:311.
PubMed  PubMed Central  Article  Google Scholar 

19.
Delebecq G, Schmidt S, Ehrhold A, Latimier M, Siano R. Revival of ancient marine dinoflagellates using molecular biostimulation. J Phycol. 2020;56:1077–89.
CAS  PubMed  Article  Google Scholar 

20.
Ribeiro S, Berge T, Lundholm N, Ellegaard M. Hundred years of environmental change and phytoplankton ecophysiological variability archived in coastal sediments. PLoS ONE. 2013;8:e61184.
CAS  PubMed  PubMed Central  Article  Google Scholar 

21.
Klouch KZ, Schmidt S, Andrieux Loyer F, Le Gac M, Hervio-Heath D, Qui-Minet ZN, et al. Historical records from dated sediment cores reveal the multidecadal dynamic of the toxic dinoflagellate Alexandrium minutum in the Bay of Brest (France). FEMS Microbiol Ecol. 2016;92:1–16.
Article  CAS  Google Scholar 

22.
Lundholm N, Ribeiro S, Godhe A, Rostgaard Nielsen L, Ellegaard M. Exploring the impact of multidecadal environmental changes on the population genetic structure of a marine primary producer. Ecol Evol. 2017;7:3132–42.
PubMed  PubMed Central  Article  Google Scholar 

23.
Moore CM, Mills MM, Arrigo KR, Berman-Frank I, Bopp L, Boyd PW, et al. Processes and patterns of oceanic nutrient limitation. Nat Geosci. 2013;6:701–10.
CAS  Article  Google Scholar 

24.
Labry C, Herbland A, Delmas D. The role of phosphorus on planktonic production of the Gironde plume waters in the Bay of Biscay. J Plankt Res. 2002;24:97–117.
CAS  Article  Google Scholar 

25.
Girault M, Arakawa H, Hashihama F. Phosphorus stress of microphytoplankton community in the western subtropical North Pacific. J Plankt Res. 2013;35:146–57.
CAS  Article  Google Scholar 

26.
Ramos JBE, Schulz KG, Voss M, Narciso Á, Müller MN, Reis FV, et al. Nutrient-specific responses of a phytoplankton community: a case study of the North Atlantic Gyre. Azores J Plankt Res. 2017;39:744–61.
Article  CAS  Google Scholar 

27.
Lin S, Litaker RW, Sunda WG. Phosphorus physiological ecology and molecular mechanisms in marine phytoplankton. J Phycol. 2016;52:10–36.
CAS  PubMed  Article  Google Scholar 

28.
Lomas MW, Swain A, Shelton R, Ammerman JW. Taxonomic variability of phosphorus stress in Sargasso Sea phytoplankton. Limnol Oceanogr. 2004;49:2303–10.
Article  Google Scholar 

29.
Wang D, Huang B, Liu X, Liu G, Wang H. Seasonal variations of phytoplankton phosphorus stress in the Yellow Sea Cold Water Mass. Acta Oceano Sin. 2014;33:124–35.
Article  CAS  Google Scholar 

30.
Cembella AD, Antia NJ, Harrison PJ. The utilization of inorganic and organic phosphorous compounds as nutrients by eukaryotic microalgae: a multidisciplinary perspective: part I. CRC Crit Rev Microbiol. 1984;10:317–91.
CAS  Article  Google Scholar 

31.
Cooper A, Bowen ID, Lloyd D. The properties and subcellular localization of acid phosphatases in the colourless alga Polytomella caeca. J Cell Sci. 1974;15:605–18.
CAS  PubMed  Google Scholar 

32.
Duhamel S, Björkman KM, Van Wambeke F, Moutin T, Karl DM. Characterization of alkaline phosphatase activity in the North and South Pacific Subtropical Gyres: Implications for phosphorus cycling. Limnol Oceanogr. 2011;56:1244–54.
CAS  Article  Google Scholar 

33.
Kang W, Wang ZH, Liu L, Guo X. Alkaline phosphatase activity in the phosphorus-limited southern Chinese coastal waters. J Environ Sci. 2019;86:38–49.
Article  Google Scholar 

34.
Girault M, Beneyton T, Pekin D, Buisson L, Bichon S, Charbonnier C, et al. High-content screening of plankton alkaline phosphatase activity in microfluidics. Anal Chem. 2018;90:4174–81. https://doi.org/10.1021/acs.analchem.8b00234.
CAS  Article  PubMed  Google Scholar 

35.
Anderson RA, Berges RA, Harrison PJ, Watanabe MM. Appendix A – recipes for freshwater and seawater media; enriched natural seawater media. In Andersen RA, editor. Algal culturing techniques. San Diego, USA: Academic; 2005. p. 429–538.

36.
Guillard RL, Ryther JH. Studies of marine planktonic diatoms. I. Cyclotella nana Hustedt, and Detonula confervacea (cleve) Gran. Can J Microbiol. 1962;8:229–39.
CAS  PubMed  Article  Google Scholar 

37.
Duffy DC, McDonald JC, Schueller OJ, Whitesides GM. Rapid prototyping of microfluidic systems in poly(dimethylsiloxane). Anal Chem. 1998;70:4974–84.
CAS  PubMed  Article  Google Scholar 

38.
Girault M, Hattori A, Kim H, Arakawa H, Matsuura K, Odaka M, et al. An on-chip imaging droplet-sorting system: a real-time shape recognition method to screen target cells in droplets with single cell resolution. Sci Rep. 2017;7:40072. https://doi.org/10.1038/srep40072.
CAS  Article  PubMed  PubMed Central  Google Scholar 

39.
Girault M, Odaka M, Kim H, Matsuura K, Terazono H, Yasuda K. Particle recognition in microfluidic applications using a template matching algorithm. JPN J Appl Phys. 2016;55. https://doi.org/10.7567/JJAP.55.06GN05.

40.
Urvoy M, Labry C, Delmas D, Creac’h L, L’Helguen S. Microbial enzymatic assays in environmental water samples: impact of Inner Filter Effect and substrate concentrations. Limnol Oceanogr Methods. 2020;18:728–38.
Article  CAS  Google Scholar 

41.
Huang Z, Terpetschnig E, You W, Haugland RP. 2-(2′-phosphoryloxyphenyl)-4-(3H)-quinazolinone derivatives as fluorogenic precipitating substrates of phosphatases. Anal Biochem. 1992;207:32–39.
CAS  PubMed  Article  Google Scholar 

42.
Girault M, Hattori A, Kim H, Matsuura K, Odaka M, Terazono H et al. Algorithm for the precise detection of single and cluster cells in microfluidic applications. Cytom Part A. 2016. https://doi.org/10.1002/cyto.a.22825.

43.
Murphy J, Riley JP. A modified single solution method for the determination of phosphate in natural waters. Anal Chim Acta. 1962;27:31–36.
CAS  Article  Google Scholar 

44.
Hoppe HG. Phosphatase activity in the sea. Hydrobiologia. 2003;493:187–200.
CAS  Article  Google Scholar 

45.
Golda-VanEeckhoutte RL, Roof LT, Needoba JA, Peterson DT. Determination of intracellular pH in phytoplankton using the fluorescent probe, SNARF, with detection by fluorescence spectroscopy. J Microbiol Methods. 2018;152:109–18.
CAS  PubMed  Article  Google Scholar 

46.
Kruskopf MM, Du Plessis S. Induction of both acid and alkaline phosphatase activity in two green-algae (chlorophyceae) in low N and P concentrations. Hydrobiologia. 2004;513:59–70.
CAS  Article  Google Scholar 

47.
Štrojsová A, Vrba J, Nedoma J, Komárková J, Znachor P. Seasonal study of extracellular phosphatase expression in the phytoplankton of a eutrophic reservoir. Eur J Phycol. 2003;38:295–306.
Article  CAS  Google Scholar 

48.
Skelton HM, Parrow MW, Burkholder JM. Phosphatase activity in the heterotrophic dinoflagellate Pfiesteria shumwayae. Harmful Algae 2006;5:395–406.
CAS  Article  Google Scholar 

49.
Nedoma J, Štrojsová A, Vrba J, Komárková J, Simek K. Extracellular phosphatase activity of natural plankton studied with ELF97 phosphate: fluorescence quantification and labelling kinetics. Environ Microbiol. 2003;5:462–72.
CAS  PubMed  Article  Google Scholar 

50.
Young EB, Tucker RC, Pansch LA. Alkaline phosphatase in freshwater Cladophora-epiphyte assemblages: regulation in response to phosphorus supply and localization. J Phycol. 2010;46:93–101.
CAS  Article  Google Scholar 

51.
Díaz-de-Quijano D, Felip M. A comparative study of fluorescence-labelled enzyme activity methods for assaying phosphatase activity in phytoplankton. A possible bias in the enzymatic pathway estimations. J Micro Meth. 2011;86:104–7.
Article  CAS  Google Scholar 

52.
Ou L, Huang B, Lin L, Hong H, Zhang F, Chen Z. Phosphorus stress of phytoplankton in the Taiwan Strait determined by bulk and single-cell alkaline phosphatase activity assays. Mar Ecol Prog Ser. 2006;327:95–106.
CAS  Article  Google Scholar 

53.
Huang B, Ou L, Wang X, Huo W, Li R, Hong H, et al. Alkaline phosphatase activity of phytoplankton in East China Sea coastal waters with frequent harmful algal bloom occurrences. Aquat Micro Ecol. 2007;49:195–206.
Article  Google Scholar 

54.
Ivančić I, Pfannkuchen M, Godrijan J, Djakovac T, Pfannkuchen DM, Korlević M, et al. Alkaline phosphatase activity related to phosphorus stress of microphytoplankton in different trophic conditions. Prog Oceanogr. 2016;146:175–86.
Article  Google Scholar 

55.
González-Gil S, Keafer B, Jovine JMR, Aguileral A, Lu S, Anderson DM. Detection and quantification of alkaline phosphatase in single cells of phosphorus-starved marine phytoplankton. Mar Ecol Prog Ser. 1998;164:21–35.
Article  Google Scholar 

56.
Dyhrman ST, Ruttenberg KC. Presence and regulation of alkaline phosphatase activity in eukaryotic phytoplankton from the coastal ocean: Implications for dissolved organic phosphorus remineralization. Limnol Oceanogr. 2006;51. https://doi.org/10.4319/lo.2006.51.3.1381.

57.
Flynn K, Jones KJ, Flynn KJ. Comparisons among species of Alexandrium (Dinophyceae) grown in nitrogen- or phosphorus-limiting batch culture. Mar Biol. 1996;126:9–18.
CAS  Article  Google Scholar 

58.
Perry MJ. Alkaline phosphatase activity in subtropical Central North Pacific waters using a sensitive fluorometric method. Mar Biol. 1972;15:113–9.
CAS  Article  Google Scholar 

59.
Dyhrman ST, Palenik B. Phosphate stress in cultures and field populations of the dinoflagellate prorocentrum minimum detected by a single-cell alkaline phosphatase assay. Appl Environ Microbiol. 1999;65:3205–12.
CAS  PubMed  PubMed Central  Article  Google Scholar 

60.
Mulholland MR, Floge S, Carpenter EJ, Capone DG. Phosphorus dynamics in cultures and natural populations of Trichodesmium spp. Mar Ecol Prog Ser. 2002;239:45–55.
CAS  Article  Google Scholar 

61.
Thomson B, Wenley J, Currie K, Hepburn C, Herndl GJ, Baltar F. Resolving the paradox: Continuous cell-free alkaline phosphatase activity despite high phosphate concentrations. Mar Chem. 2019;214:103671.
CAS  Article  Google Scholar 

62.
Foster RA, Sztejrenszus S, Kuypers MMM. Measuring carbon and N2 fixation in field populations of colonial and free-living unicellular cyanobacteria using nanometer-scale secondary ion mass spectrometry. J Phycol. 2013;49:502–16.
CAS  PubMed  Article  Google Scholar 

63.
Dyhrman ST, Palenik B. Characterization of ectoenzyme activity and phosphate-regulated proteins in the coccolithophorid Emiliania huxleyi. J Plankton Res. 2003;25:1215–25.
CAS  Article  Google Scholar 

64.
Oh SJ, Yamammoto T, Kataoka Y, Matsuda O, Matsuyama Y, Katani Y. Utilization of dissolved organic phosphorus by the two toxic dinoflagellates, Alexandrium tamarense and Gymnodinium catenatum (Dinophyceae). Fish Sci. 2002;68:416–24.
CAS  Article  Google Scholar 

65.
Jauzein C, Labry C, Youenou A, Quéré J, Delmas D, Collos Y. Growth and phosphorus uptake by the toxic dinoflagellate Alexandrium catenella (Dinophycea) in response to phosphate limitation. J Phycol. 2010;46:926–36.
CAS  Article  Google Scholar 

66.
Elgavish A, Halmann M, Berman T. A comparative study of phosphorus utilization and storage in batch cultures of Peridinium cinctum, Pediastrum duplex and Cosmarium sp., from Lake Kinneret (Israel). Phycologia. 1982;21:47–54.
CAS  Article  Google Scholar 

67.
Flynn K, Franco JM, Fernandez P, Reguera B, Zapata M, Wood G, et al. Changes in toxin content, biomass and pigments of the dinoflagellate Alexandrium minutum during nitrogen refeeding and growth into nitrogen or phosphorus stress. Mar Ecol Prog Ser. 1994;111:99–109.
CAS  Article  Google Scholar 

68.
Ou L, Wang D, Huang B, Hong H, Qi Y, Lu S. Comparative study of phosphorus strategies of three typical harmful algae in Chinese coastal waters. J Plankton Res. 2008;30:1007–17.
CAS  Article  Google Scholar 

69.
Droop MR. The nutrient status of algal cells in continuous culture. J Mar Biol Ass UK. 1974;54:825–55.
CAS  Article  Google Scholar 

70.
Bechemin C, Grzebyk D, Hachame F, Hummert C, Maestrini S. Effect of different nitrogen/phosphorus nutrient ratios on the toxin content in Alexandrium minutum. Aquat Micro Ecol. 1990;20:157–65.
Article  Google Scholar 

71.
Labry C, Erard–Le Denn E, Chapelle A, Fauchot J, Youenou A, Crassous MP, et al. Competition for phosphorus between two dinoflagellates: A toxic Alexandrium minutum and a non-toxic Heterocapsa triquetra. J Exp Mar Biol Ecol. 2008;358:124–35.
CAS  Article  Google Scholar 

72.
Chapelle A, Labry C, Sourisseau M, Lebreton C, Youenou A, Crassous MP. Alexandrium minutum growth controlled by phosphorus An applied model. J Mar Syst. 2010;83:181–91.
Article  Google Scholar 

73.
Sakshaug E, Granéli E, Elbrächter M, Kayser H. Chemical composition and alkaline phosphatase activity of nutrient-saturated and P-deficient cells of four marine dinoflagellates. J Exp Mar Biol Ecol. 1984;11:241–54.
Article  Google Scholar 

74.
Lirdwitayaprasit T, Okaichi T, Montani S, Ochi T, Anderson DM. Changes in cell chemical con~position during the life cycle of Scrippsiella trochoidea (Dinophyceae). J Phycol. 1990;26:299–306.
CAS  Article  Google Scholar 

75.
Qi H, Wang J, Wang Z. A comparative study of maximal quantum yield of photosystem II to determine nitrogen and phosphorus limitation on two marine algae. J Sea Res. 2013;80:1–11.
Article  Google Scholar 

76.
Simon N, Cras AL, Foulon E, Lemée R. Diversity and evolution of marine phytoplankton. C R Biol. 2009;332:159–70.
PubMed  Article  Google Scholar 

77.
Rengefors K, Kremp A, Reusch TBH, Wood AM. Genetic diversity and evolution in eukaryotic phytoplankton: revelations from population genetic studies. J Plankton Res. 2017;39:165–79.
Google Scholar 

78.
Bendif EM, Nevado B, Wong ELY, Wong EL, Hagino K, Probert I, et al. Repeated species radiations in the recent evolution of the key marine phytoplankton lineage Gephyrocapsa. Nat Commun. 2019;10:4234.
PubMed  PubMed Central  Article  CAS  Google Scholar 

79.
Thornton DCO. Individuals clones or groups? Phytoplankton behaviour and units of selection. Ethol Ecol Evol. 2002;14:165–73.
Article  Google Scholar 

80.
Gerecht A, Romano G, Lanora A, d’Ippolito G, Cutignano A, Fontana A. Plasticity of Oxylipin metabolism among clones of the marine diatom Skeletonema marinoi (Bacillariophyceae). J Phycol. 2011;47:1050–6.
CAS  PubMed  Article  Google Scholar 

81.
Lim PT, Leaw CP, Usup G, Kobiyama A, Koike K, Ogata T. Effects of light and temperature on growth, nitrate uptake, and toxin production of two tropical dinoflagellates: Alexandrium tamiyavanichii and Alexandrium minutum (Dinophyceae). J Phycol. 2006;42:786–99.
CAS  Article  Google Scholar 

82.
Van Mooy BA, Fredricks HF, Pedler BE, Dyhrman ST, Karl DM, Koblížek M, et al. Phytoplankton in the ocean use non-phosphorus lipids in response to phosphorus scarcity. Nature. 2009;458:69–72.
PubMed  Article  CAS  Google Scholar 

83.
Galbraith AD, Martiny AC. Simple mechanism for marine nutrient stoichiometry. Proc Natl Acad Sci USA. 2015;112:8199–204.
CAS  PubMed  Article  Google Scholar 

84.
Berge T, Daugbjerg N, Hansen PJ. Isolation and cultivation of microalgae select for low growth rate and tolerance to high pH. Harmful Algae. 2012;20:101–10.
CAS  Article  Google Scholar  More

Composition of chemical constituents and bacterial communities in AGS and feces indicates separate, unique odor profiles
The gas chromatography–mass spectrometry analyses revealed that AGS volatiles of wild and captive pandas were comprised of a multicomponent blend of 30–50 chemical compounds, including fatty acids, aldehydes, ketones, aliphatic ethers, amides, aromatics, alcohols, steroids and squalene (Fig. 2a and Supplementary Table S2). These compounds are typical components of chemosignals across species due to their volatility, detectability and other characteristics facilitating chemoreception [3, 26, 32]. By contrast, feces contained mostly fatty acid ethyl ester, and a small number and quantity of fatty acids, amides, steroids and indole (Fig. 2b and Supplementary Table S3). Our results show that the relative abundance of steroids, aldehydes and fatty acids were remarkably higher in AGS than in feces (Fig. 3a), and the number of chemical components of aldehydes, fatty acids, and ketones in AGS was also significantly higher than found in feces (Fig. 3b). These results indicate that the chemical constituents of AGS are much better suited for chemosignaling than those from feces.
Fig. 2: Representative ion chromatograms of samples in giant pandas.

a Anogenital gland secretions (AGS). b Feces.

Full size image

Fig. 3: Differences in chemical compounds of anogenital gland secretions (AGS) and feces in giant pandas, and the differences in microbial communities, KEGG and contribution bacteria for lipid metabolism.

a A heat map of the mean relative abundance of the chemical compounds. b A heat map of the number chemical components. Differences in the microbial communities as a function of providence (captive/wild) and source (feces/AGS) at the c phylum and d genus level. e PCoA clustering results of samples from different groups. f Hierarchical clustering analysis of the samples, clearly indicating two branches for AGS and fecal samples. g Six differentially represented pathways in lipid metabolism and the Linear discriminant analysis (LDA) score. h Prevalence of enzymes involved in lipid metabolism as a function of phylum and family in AGS of giant pandas. i The contribution of different bacteria at genus level to lipid metabolism. WPF: wild panda feces, CPF: captive panda feces, WPAG: wild panda AG, CPAG: captive panda AG.

Full size image

The composition of bacterial communities in AGS and feces was markedly different at the phylum (Fig. 3c) and genus levels (Fig. 3d), based on taxonomic classifications of predicted gene sequences. Principal Co-ordinates Analysis (PCoA) (Fig. 3e) and hierarchical clustering analyses (Fig. 3f) revealed cluster patterns based on provenance (captive/wild) and sample type (AGS/fecal). Notably, the microbiota composition of AGS from different individuals or living environments was more similar than were AGS and fecal samples from the same individuals. Actinobacteria (X2 = 26.33, P  More