Ecology

Different sampling strategies yield different microbial communities
The sampling strategies compared in this study (homogenizing tissue before subsampling and homogenizing tissue after subsampling) are common methods found in the literature for characterizing plant-associated microbial communities14,23,26,29. Procrustes analyses and community overlap between sampling strategies demonstrated that different strategies can capture disparate microbial communities within plants, with the extent of these differences depending on the community targeted and plant tissue type sampled. In FFE as well as bacterial and non-AM fungal communities in roots, subsamples from the same plant resulted in completely different sets of species recovered, illustrating the severe undersampling that is inherent to each of these strategies. With these sampling strategies, we are undoubtedly sacrificing power and accuracy to characterize the subtler aspects of plant microbiome interactions, despite often seeing community differences across landscapes, treatments or seasons.
Richness was higher when homogenizing before subsampling for bacteria only, despite differences observed in composition for all groups. It is perhaps surprising that homogenizing plant tissues before subsampling did not recover more species than homogenizing after subsampling for fungi as well, because with the former approach, more plant tissue is initially represented. Indeed, a previous study showed that sample pooling or homogenizing before subsampling resulted in a higher richness of soil fungi compared to equally sized individual samples50. In Song et al. (2015)50 they also found that multiple individual subsamples, rather than the single homogenized subsample, resulted in higher richness. This may suggest that the scale at which we are physically able to break down the particle size of plant tissues, as opposed to soil, is not always fine enough to sufficiently homogenize the fungi within. Because of this, plant-associated microbial communities may require a greater sampling effort than soil microbes. Additionally, the removal of low-abundant SVs did not result in differences in richness between the two sampling strategies for any microbial group, suggesting that neither strategy is better at capturing rare species. Although this study was performed only on milkweed plants, we believe that these results are applicable to other plant species as well. The richness reported here is similar to other studies of plant-associated microbes (e.g.51,52), indicating that differences in subsamples were not due to extreme richness of milkweed-associated microbes.
Microbial diversity should inform sampling effort
The higher congruency that we saw between sampling strategies for AMF compared to other microbial communities may be due to the differences in their local and global estimated richness. While the global number of AMF species has been estimated in the hundreds to low thousands42,53, global estimates of fungal species in general range in the millions54,55. A recent global estimate of bacterial richness suggests similar scales56. In this study specifically, AMF had the lowest total SV richness and the greatest similarity between sampling strategies, while foliar fungal endophytes had the highest SV richness, and the lowest overlap of SVs between strategies. Since the amount of tissue sampled was equal for all microbial communities, the sampling effort was likely much higher for AMF (relative to the whole AMF community), than it was for bacteria and non-AM fungi. Consequently, with each sample we are likely sampling a much larger proportion of true AMF species richness.
Even though the estimated total community richness was highest for foliar fungi, the average estimated richness per individual plant was highest for AMF. This suggests that similar AMF SVs re-occurred across all plants with low species turnover. On the other hand, fungi in leaves had lower average richness per plant (Fig. 4, Supplementary Fig. S3), but the highest total richness, meaning that there was higher turnover of FFE species among plants sampled. These results may be a direct reflection of the overall community richness of the different microbial groups as well as their ability to spread and co-occur within plants. Based on these patterns, more individual plants and a greater sampling effort within individuals are likely needed to characterize FFE communities compared to AMF communities.
Rare SVs contribute to variation among subsamples
Our results show that low abundant, rare SVs largely contributed to the differences seen between sampling strategies. Even AMF communities, which were already similar, increased in overlap by 50% between strategies after low abundant SVs (represented by  More

1.
Hallmann, C. A. et al. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS ONE 12, e0185809 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 
2.
Sánchez-Bayo, F. & Wyckhuys, K. A. G. Worldwide decline of the entomofauna: a review of its drivers. Biol. Conserv. 232, 8–27 (2019).
Article  Google Scholar 

3.
van Klink, R. et al. Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances. Science 368, 417–420 (2020).
ADS  PubMed  Article  CAS  PubMed Central  Google Scholar 

4.
Losey, E. J. & Vaughan, M. The economic value of ecological services provided by insects. Bioscience 56, 311 (2006).
Article  Google Scholar 

5.
Bar-On, Y. M., Phillips, R. & Milo, R. The biomass distribution on Earth. Proc. Natl. Acad. Sci. USA 115, 6506–6511 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

6.
Yang, L. H. & Gratton, C. Insects as drivers of ecosystem processes. Curr. Opin. Insect Sci. 2, 26–32 (2014).
PubMed  Article  PubMed Central  Google Scholar 

7.
Hunter, M. D. Insect population dynamics meets ecosystem ecology: effects of herbivory on soil nutrient dynamics. Agric. For. Entomol. 3, 77–84 (2001).
Article  Google Scholar 

8.
Nichols, E. et al. Ecological functions and ecosystem services provided by Scarabaeinae dung beetles. Biol. Conserv. 141, 1461–1474 (2008).
Article  Google Scholar 

9.
Lobry de Bruyn, L. A. & Conacher, A. J. The role of termites and ants in soil modification: a review. Aust. J. Soil Res. 28, 55–93 (1990).
Google Scholar 

10.
Shukla, R. K., Singh, H., Rastogi, N. & Agarwal, V. M. Impact of abundant Pheidole ant species on soil nutrients in relation to the food biology of the species. Appl. Soil Ecol. 71, 15–23 (2013).
Article  Google Scholar 

11.
López-Hernández, D. Nutrient dynamics (C, N and P) in termite mounds of Nasutitermes ephratae from savannas of the Orinoco Llanos (Venezuela). Soil Biol. Biochem. 33, 747–753 (2001).
Article  Google Scholar 

12.
Nkem, J. N., Lobry de Bruyn, L. A., Grant, C. D. & Hulugalle, N. R. The impact of ant bioturbation and foraging activities on surrounding soil properties. Pedobiologia 44, 609–621 (2000).
Article  Google Scholar 

13.
Jouquet, P., Traoré, S., Choosai, C., Hartmann, C. & Bignell, D. Influence of termites on ecosystem functioning. Ecosystem services provided by termites. Eur. J. Soil Biol. 47, 215–222 (2011).
Article  Google Scholar 

14.
Frost, C. J. & Hunter, M. D. Insect canopy herbivory and frass deposition affect soil nutrient dynamics and export in oak mesocosms. Ecology 85, 3335–3347 (2004).
Article  Google Scholar 

15.
Calderon-Cortes, N. et al. Ecosystem engineering and manipulation of host plant tissues by the insect borer Oncideres albomarginata chamela. J. Insect Physiol. 84, 128–136 (2016).
CAS  PubMed  Article  Google Scholar 

16.
Barton, P. S., Cunningham, S. A., Lindenmayer, D. B. & Manning, A. D. The role of carrion in maintaining biodiversity and ecological processes in terrestrial ecosystems. Oecologia 171, 761–772 (2013).
ADS  PubMed  Article  Google Scholar 

17.
Danell, K., Aux, D. B. & Bráthen, K. A. Effect of muskox carcasses on nitrogen concentration in tundra vegetation. Artic 55, 389–392 (2002).
Google Scholar 

18.
Melis, C. et al. Soil and vegetation nutrient response to bison carcasses in Bialowieza primeval forest Poland. Ecol. Res. 22, 807–813 (2007).
ADS  CAS  Article  Google Scholar 

19.
Bump, J. K., Peterson, R. O. & Vucetich, J. A. Wolves modulate soil nutrient heterogeneity and foliar nitrogen by configuring the distribution of ungulate carcasses. Ecology 90, 3159–3167 (2009).
PubMed  Article  Google Scholar 

20.
Benninger, L. A., Carter, D. O. & Forbes, S. L. The biochemical alteration of soil beneath a decomposing carcass. Forensic Sci. Int. 180, 70–75 (2008).
CAS  PubMed  Article  Google Scholar 

21.
Anderson, G. S. & VanLaerhoven, S. L. Initial studies on insect succession on carrion in Southwestern British Columbia. J. Forensic Sci. 41, 617–625 (1996).
Article  Google Scholar 

22.
Matuszewski, S., Bajerlein, D., Konwerski, S. & Szpila, K. Insect succession and carrion decomposition in selected forests of Central Europe. Part 3: succession of carrion fauna. Forensic Sci. Int. 207, 150–163 (2011).
PubMed  Article  Google Scholar 

23.
Barton, P. S., Evans, M. J., Pechal, J. L. & Benbow, M. E. Necrophilous insect dynamics at small vertebrate carrion in a temperate eucalypt woodland. J. Med. Entomol. 54, 964–973 (2017).
PubMed  Article  Google Scholar 

24.
Benbow, M. E., Lewis, A. J., Tomberlin, J. K. & Pechal, J. L. Seasonal necrophagous insect community assembly during vertebrate carrion decomposition. J. Med. Entomol. 50, 440–450 (2013).
CAS  PubMed  Article  Google Scholar 

25.
Payne, J. A. A summer carrion study of the baby pig Sus Scrofa Linnaeus. Ecology 46, 592–602 (1965).
Article  Google Scholar 

26.
Kočárek, P. Decomposition and Coleoptera succession on exposed carrion of small mammal in Opava, the Czech Republic. Eur. J. Soil Biol. 39, 31–45 (2003).
Article  Google Scholar 

27.
Bornemissza, G. F. An analysis of arthropod succession in carrion and the effects of its decomposition on the soil fauna. J. Zool. 5, 1–12 (1957).
Google Scholar 

28.
Parmenter, R. R. & MacMahon, J. A. Carrion decomposition and nutrient cycling in a semiarid shrub–steppe ecosystem. Ecol. Monogr. 79, 637–661 (2009).
Article  Google Scholar 

29.
Pechal, J. L., Benbow, M. E., Crippen, T. L., Tarone, A. M. & Tomberlin, J. K. Delayed insect access alters carrion decomposition and necrophagous insect community assembly. Ecosphere 5, 1–21 (2014).
Article  Google Scholar 

30.
Pukowski, E. Ökologische untersuchungen an necrophorus. Z. Morphol. Oekol. 27, 518–586 (1933).
Article  Google Scholar 

31.
Scott, M. P. The ecology and behavior of burying beetles. Annu. Rev. Entomol. 43, 595–618 (1998).
CAS  PubMed  Article  Google Scholar 

32.
Milne, L. J. & Milne, M. The social behavior of burying beetles. Sci. Am. 235, 84–89 (1976).
Article  Google Scholar 

33.
Müller, J. K. & Eggert, A. K. Paternity assurance by ‘helpful’ males: adaptations to sperm competition in burying beetles. Behav. Ecol. Sociobiol. 24, 245–249 (1989).
Article  Google Scholar 

34.
House, C. M. et al. The evolution of repeated mating in the burying beetle Nicorphorus vespilloides. Evolution 62, 2004–2014 (2008).
PubMed  Article  Google Scholar 

35.
Pettinger, A. M., Steiger, S., Müller, J. K., Sakaluk, S. K. & Eggert, A. K. Dominance status and carcass availability affect the outcome of sperm competition in burying beetles. Behav. Ecol. 22, 1079–1087 (2011).
Article  Google Scholar 

36.
Rozen, D. E., Engelmoer, D. J. P. & Smiseth, P. T. Antimicrobial strategies in burying beetles breeding on carrion. Proc. Natl. Acad. Sci. USA 105, 17890–17895 (2008).
ADS  CAS  PubMed  Article  Google Scholar 

37.
Arce, A. N., Johnston, P. R., Smiseth, P. T. & Rozen, D. E. Mechanisms and fitness effects of antibacterial defences in a carrion beetle. J. Evol. Biol. 25, 930–937 (2012).
CAS  PubMed  Article  Google Scholar 

38.
Hall, C. L. et al. Inhibition of microorganisms on a carrion breeding resource: the antimicrobial peptide activity of burying beetle (Coleoptera: Silphidae) oral and anal secretions. Environ. Entomol. 40, 669–678 (2011).
PubMed  Article  Google Scholar 

39.
Duarte, A., Welch, M., Swannack, C., Wagner, J. & Kilner, R. M. Strategies for managing rival bacterial communities: lessons from burying beetles. J. Anim. Ecol. 87, 414–427 (2018).
PubMed  Article  Google Scholar 

40.
Trumbo, S. T. Reproductive benefits and the duration of paternal care in a biparental burying beetle Necrophorus orbicollis. Behaviour 117, 82–105 (1991).
Article  Google Scholar 

41.
Trumbo, S. T. Fate of mouse carcasses in a Northern Woodland. Ecol. Entomol. 41, 737–740 (2016).
Article  Google Scholar 

42.
Hoback, W. W., Freeman, L., Payton, M. & Peterson, B. C. Burying beetle (Coleoptera: Silphidae: Nicrophorus fabricius) brooding improves soil fertility. Coleopt. Bull. 74, 427–433 (2020).
Article  Google Scholar 

43.
Benbow, M. E. et al. Necrobiome framework for bridging decomposition ecology of autotrophically and heterotrophically derived organic matter. Ecol. Monogr. 89, 1–26 (2018).
Google Scholar 

44.
Carter, D. O., Yellowlees, D. & Tibbett, M. Cadaver decomposition in terrestrial ecosystems. Naturwissenschaften 94, 12–24 (2007).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

45.
Carter, D. O., Yellowlees, D. & Tibbett, M. Temperature affects microbial decomposition of cadavers (Rattus rattus) in contrasting soils. Appl. Soil Ecol. 40, 129–137 (2008).
Article  Google Scholar 

46.
Meyer, J., Anderson, B. & Carter, D. O. Seasonal variation of carcass decomposition and gravesoil chemistry in a cold (Dfa) climate. J. Forensic Sci. 58, 1175–1182 (2013).
PubMed  Article  PubMed Central  Google Scholar 

47.
Metcalf, J. L. et al. Microbial community assembly and metabolic function during mammalian corpse decomposition. Science 351, 158–162 (2016).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

48.
Chen, Y. A., Forschler, B. T. & Nielsen, U. Elemental concentrations in the frass of saproxylic insects suggest a role in micronutrient cycling. Ecosphere 7, 1–13 (2016).
Google Scholar 

49.
Couture, J. J. & Lindroth, R. L. Atmospheric change alters frass quality of forest canopy herbivores. Arthropod. Plant. Interact. 8, 33–47 (2014).
Article  Google Scholar 

50.
Metcalf, J. L. et al. A microbial clock provides an accurate estimate of the postmortem interval in a mouse model system. Elife 2013, 1–19 (2013).
Google Scholar 

51.
Xu, S., Liu, L. L. & Sayer, E. J. Variability of above-ground litter inputs alters soil physicochemical and biological processes: a meta-analysis of litterfall-manipulation experiments. Biogeosciences 10, 7423–7433 (2013).
ADS  Article  CAS  Google Scholar 

52.
Wilson, W. A. et al. Regulation of glycogen metabolism in yeast and bacteria. FEMS Microbiol. Rev. 34, 952–985 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

53.
Achbergerová, L. & Nahálka, J. Polyphosphate—an ancient energy source and active metabolic regulator. Microb. Cell Fact. 10, 1–14 (2011).
Article  CAS  Google Scholar 

54.
Kornberg, A. Inorganic polyphosphate: toward making a forgotten polymer unforgettable. J. Bacteriol. 177, 491–496 (1995).
CAS  PubMed  PubMed Central  Article  Google Scholar 

55.
Wilkinson, J. F. Carbon and energy storage in bacteria. J. Gen. Microbiol. 32, 171–176 (1963).
CAS  PubMed  Article  Google Scholar 

56.
Paul, E. A. Soil Microbiology, Ecology and Biochemistry (Academic Press, Cambridge, 2014).
Google Scholar 

57.
DeVault, T. L., Brisbin, I. L. Jr. & Rhodes, O. E. Jr. Factors influencing the acquisition of rodent carrion by vertebrate scavengers and decomposers. Can. J. Zool. 82, 502–509 (2004).
Article  Google Scholar 

58.
DeVault, T. L., Rhodes, O. E. Jr. & Shivik, J. A. Scavenging by vertebrates: behavioral, ecological and evolutionary perspectives on an important energy transfer pathway in terrestrial ecosystems. Oikos 102, 225–234 (2003).
Article  Google Scholar 

59.
DeVault, T. L. & Rhodes, O. E. Identification of vertebrate scavengers of small mammal carcasses in a forested landscape. Acta Theriol. 47, 185–192 (2002).
Article  Google Scholar 

60.
Smith, J. B., Laatsch, L. J. & Beasley, J. C. Spatial complexity of carcass location influences vertebrate scavenger efficiency and species composition. Sci. Rep. 7, 1–8 (2017).
Article  CAS  Google Scholar 

61.
Wilson, D. S. & Fudge, J. Burying beetles: intraspecific interactions and reproductive success in the field. Ecol. Entomol. 9, 195–203 (1984).
Article  Google Scholar 

62.
Keller, M. L., Howard, D. R. & Hall, C. L. Spatiotemporal niche partitioning in a specious silphid community (Coleoptera: Silphidae Nicrophorus). Sci. Nat. 106, 1–12 (2019).
ADS  Article  CAS  Google Scholar 

63.
Owings, C. G. & Picard, C. J. Temporal survey of a carrion beetle (Coleoptera:Silphidae) community in Indiana. Proc. Indiana Acad. Sci. 124, 124–128 (2015).
Google Scholar 

64.
Snyder, D. P. Survival rates, longevity, and population fluctuations in the white-footed mouse, Peromyscus leucopus Southeastern Michigan. Museum Zool. Universty Michigan 95, 1–33 (1956).
Google Scholar 

65.
Stephens, R. B., Hocking, D. J., Yamasaki, M. & Rowe, R. J. Synchrony in small mammal community dynamics across a forested landscape. Ecography 40, 1198–1209 (2017).
Article  Google Scholar 

66.
Leasure, D. R., Rupe, D. M., Phillips, E. A., Opine, D. R. & Huxel, G. R. Efficient new above-ground bucket traps produce comparable data to that of standard transects for the endangered American Burying Beetle, Nicrophorus americanus Olivier (Coleoptera: Silphidae). Coleopt. Bull. 66, 209–218 (2012).
Article  Google Scholar 

67.
Bedick, J. C., Ratcliffe, B. C. & Higley, L. G. A new sampling protocol for the endangered American Burying Beetle, Nicrophorus americanus Olivier (Coleoptera: Silphidae). Coleopt. Bull. 58, 57–70 (2004).
Article  Google Scholar 

68.
Trumbo, S. T. Reproductive success, phenology and biogeography of burying beetles (Silphidae, Nicrophorus). Am. Midl. Nat. 124, 1–11 (1990).
Article  Google Scholar 

69.
Trumbo, S. T. Nesting failure in burying beetles and the origin of communal associations. Evol. Ecol. 9, 125–130 (1995).
Article  Google Scholar 

70.
Braman, R. S. & Hendrix, S. A. Nanogram nitrite and nitrate determination in environmental and biological materials by vanadium (III) reduction with chemiluminescence detection. Anal. Chem. 61, 2715–2718 (1989).
CAS  PubMed  Article  Google Scholar 

71.
Sims, G. K., Ellsworth, T. R. & Mulvaney, R. L. Microscale determination of inorganic nitrogen in water and soil extracts. Commun. Soil Sci. Plant Anal. 26, 303–316 (1995).
CAS  Article  Google Scholar 

72.
Contosta, A. R., Frey, S. D. & Cooper, A. B. Seasonal dynamics of soil respiration and N mineralization in chronically warmed and fertilized soils. Ecosphere 2, 1–21 (2011).
Article  Google Scholar 

73.
Vance, E. D., Brookes, P. C. & Jenkinson, D. S. Microbial biomass measurements in forest soils: the use of the chloroform fumigation-incubation method in strongly acid soils. Soil Biol. Biochem. 19, 697–702 (1987).
CAS  Article  Google Scholar 

74.
Vance, E. D., Brookes, P. C. & Jenkinson, D. S. An extraction method for measuring soil microbial biomass C. Soil Biol. Biochem. 19, 703–707 (1987).
CAS  Article  Google Scholar 

75.
Brookes, P. D. C., Landman, A., Pruden, G. & Jenkinson, D. S. Chloroform fumigation and the release of soil nitrogen: a rapid direct extraction method to measure microbial biomass nitrogen in soil. Soil Biol. Biochem. 17, 837–842 (1985).
CAS  Article  Google Scholar 

76.
Bligh, E. G. & Dyer, W. J. A rapid method of total lipid extraction and purification. Can. J. Biochem. Physiol. 37, 911–917 (1959).
CAS  PubMed  Article  Google Scholar 

77.
White, D. C., Davis, W. M., Nickels, J. S., King, J. D. & Bobbie, R. J. Determination of the sedimentary microbial biomass by extractible lipid phosphate. Oecologia 40, 51–62 (1979).
ADS  CAS  PubMed  Article  Google Scholar 

78.
Guckert, J. B., Antworth, C. P., Nichols, P. D. & White, D. C. Phospholipid, ester-linked fatty acid profiles as reproducible assays for changes in prokaryotic community structure of estuarine sediments. FEMS Microbiol. Lett. 31, 147–158 (1985).
CAS  Article  Google Scholar 

79.
Bardgett, R. D., Hobbs, P. J. & Frostegårde, A. Changes in soil fungal:bacterial biomass ratios following reductions in the intensity of management of an upland grassland. Biol. Fertil. Soils 22, 261–264 (1996).
Article  Google Scholar 

80.
Bååth, E. The use of neutral lipid fatty acids to indicate the physiological conditions of soil fungi. Microb. Ecol. 45, 373–383 (2003).
PubMed  Article  CAS  Google Scholar 

81.
Ekelund, F., Olsson, S. & Johansen, A. Changes in the succession and diversity of protozoan and microbial populations in soil spiked with a range of copper concentrations. Soil Biol. Biochem. 35, 1507–1516 (2003).
CAS  Article  Google Scholar 

82.
Leckie, S. E., Prescott, C. E., Grayston, S. J., Neufeld, J. D. & Mohn, W. W. Characterization of humus microbial communities in adjacent forest types that differ in nitrogen availability. Microb. Ecol. 48, 29–40 (2004).
CAS  PubMed  Article  Google Scholar 

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

84.
Oksanen, J. et al. Vegan: community ecology package (2019). http://cran.r-project.org/.

85.
Hervé, M. RVAideMemoire: Testing and plotting procedures for biostatistics version 0.9-75 from CRAN. R package version 0.9-75 (2020).

86.
Fox, J. & Weisberg, S. An R Companion to Applied Regression (SAGE Publications Inc, Thousand Oaks, 2018).
Google Scholar  More

In this section, the heterogeneity of the landscape is introduced by assuming that its profile can be written as (Psi (x) = a + psi (x)), where (psi (x)) represents the spatial variations of the environment around a reference level a.
The results that we will present were obtained through theoretical and numerical techniques. The theoretical approach is based on the mode linear stability analysis discussed in the previous section. Numerical integration of Eq. (3), starting from a homogeneous state plus a random perturbation, was performed following an explicit forward-time-centered-space scheme, with boundary conditions suitably chosen for each case (see Supplementary Information for details).
Refuge
As a paradigm of a heterogeneous environment with sharp borders, we first consider that the spatial variations around the reference level a are given by

$$begin{aligned} psi (x) = – A[1- Theta (L/2 -|x|)] ,, end{aligned}$$
(6)

where (Theta) is the Heaviside step function. With (A >0), it represents a refuge of size L with growth rate a immersed in a less viable environment with growth rate (a-A). In a laboratory situation, this can be constructed by means of a mask delimiting a region that protects organisms from some harmful agent, for instance, shielding bacteria from UV radiation26. In natural environments, this type of localized disturbance appears due to changes in the geographical and local climate conditions27, or even engineered by other species38.
In Homogeneous landscapes section, we have seen that the uniform distribution is intrinsically stable when (q ge 0). In contrast, when there are heterogeneities in (Psi (x)), spatial structures can emerge even if (qge 0), as illustrated in Fig. 3 for the case (D=0.01).
Figure 3

Stationary population density (rho _s) vs. x in a refuge. This heterogeneous environment is defined by Eq. (6), with (a=b=1), (A=10^{-3}) and (L=10). The vertical lines indicate the refuge boundaries. We used the kernel (gamma _q(x)), with (q=0.1) and (ell =2), and two different values of D. Symbols are results from numerical integration of Eq. (3) under periodic boundary conditions, and solid lines from the small-A approximation given by Eq. (8), in excellent agreement with the exact numerical solution. Recall that, in a homogeneous environment, no oscillations appear for (q ge 0).

Full size image

In the limit of weak heterogeneity, i.e., under the condition (|psi (x)|/a ll 1), we obtain an approximate analytical solution assuming that the steady solution of Eq. (3) can be expressed in terms of a small deviation (varepsilon _s(x)) around the homogeneous state (rho _0=a/b). Then, we substitute (rho _s(x)=rho _0+varepsilon _s(x)) into the stationary form of Eq. (3), discard terms of order equal or higher than (mathcal{O}(varepsilon ^2, Avarepsilon ,A^2)), and Fourier transform, obtaining

$$begin{aligned} tilde{varepsilon }_s(k) = dfrac{ rho _0 tilde{psi }(k)}{-lambda (k)},, end{aligned}$$
(7)

where (lambda (k)) was already defined in Eq. (5) and (tilde{psi }(k)) is the Fourier transform of the small fluctuations in the landscape quality, which for the case of Eq. (6) is (tilde{psi }(k)= A[2sin (Lk/2)/k -2pi delta (k)]).
Finally, assuming that (lambda (k^star )0) and (bar{x}rightarrow infty), lilac), decaying oscillations ((bar{k} >0) and finite (bar{x}), orange), and pure exponential decay ((bar{k}=0) and finite (bar{x}), gray). In (a), the dashed and dotted lines correspond to ({k_i}=0) and ({k_r}=0), respectively, where ({k_r}) and ({k_i}) are the real and imaginary parts of the zeros of (lambda(k)), with the smallest positive imaginary part. In (b)-(d), symbols correspond to measurements of numerical profiles, according to Fig. 4, and solid lines correspond to the prediction in Eq. (10) (theoretical 1). Thin dashed lines correspond to the harmonic estimate (theoretical 2) given by Eq. (14).

Full size image

To perform a theoretical prediction of (bar{k}) and (bar{x}), within the linear approximation, we consider that these oscillation parameters should be related to the poles of the integrand (mathrm{e}^{ikx} tilde{psi }(k) /[-lambda (k)]) in the expression for the inverse Fourier transform that provides the solution, according to Eq. (8). As far as the external field (psi (x)) does not introduce non-trivial poles, like in the case of a Heaviside step function ((tilde{psi }(k) sim 1/k)), only the zeros of the complex extension of (lambda (k)) matter. The dominant (more slowly decaying mode) is given by the complex poles (k=pm k_r + i k_i) ( (k_r >0)) with minimal positive imaginary part that, except for amplitude and phase constants, will approximately provide patterns of the form (mathrm{e}^{-k_i x} cos (k_r x)), allowing the identifications

$$begin{aligned} bar{k} = k_rquad text { and }quad 1/bar{x}= k_i. end{aligned}$$
(10)

This theoretical prediction42 is in very good agreement with the results of numerical simulations, as shown in Fig. 5, explaining the observed regimes.
Moreover, the modes that persist beyond the interface have relatively small amplitudes, so that the system response is approximately linear in this region.
Lastly, recall that this analysis assumes mode stability ((lambda (k)0), the system is intrinsically unstable, with the poles having null imaginary part (lying on the real axis). Nevertheless, the initially fastest growing mode, given by the maximum of (lambda (k)), tends to remain the dominant one in the long term41, yielding (bar{k} simeq k^star) for the sustained oscillations ((bar{x} rightarrow infty)).
In order to obtain further insights, it is useful to consider the response function (tilde{R}(k)) that, from Eq. (7), is

$$begin{aligned} tilde{R}(k) equiv frac{|tilde{varepsilon }_s(k)|^2}{|tilde{psi }(k)|^2} = frac{rho _0^2}{lambda ^{2}(k)} ,. end{aligned}$$
(11)

Despite missing some of the dynamical information contained in the phase of (lambda (k)), it can provide a more direct estimation of the observed parameters than through calculation of the poles. In order to perform this estimation, we resort to the response function of a driven damped linear oscillator43 described by the equation (varepsilon _H”(x)+2zeta k_0varepsilon _H'(x)+k_0^2varepsilon _H(x)= f(x)). We have

$$begin{aligned} tilde{R}_H(k) equiv frac{|tilde{varepsilon }_H(k)|^2}{|tilde{f}(k)|^2}= frac{1}{|lambda _H(k)|^{2}} = frac{1}{(k^2-k_0^2)^2+4zeta ^2k_0^2k^2}, , end{aligned}$$
(12)

with (-lambda _H(k) = -k^2 + i2zeta k_0 k + k_0^2), whose zeros (poles of (1/lambda _H(k))) are (k= pm k_r + i k_i = k_0 (pm sqrt{1-zeta ^2}+izeta )), where (k_0) is the natural mode and (zeta) is the damping coefficient. Note that, under a step forcing (f(x)=k_0^2 Theta (x)), which simulates our present setting, those poles carry the essential information of the damped-oscillation solution, given by (tilde{varepsilon }_H(k) = tilde{f}(k)/[-lambda _H(k)]), where (tilde{f}(k)=k_0^2(pi delta (k) -i/k)). In the underdamped case ((zeta 1), when the zeros of (lambda (k)) are pure imaginary with (k_i=k_0(zeta pm sqrt{zeta ^2-1})). The connection between the poles of (tilde{R}_H(k)) and the dynamic solution of the driven harmonic oscillator is possible because, as previously discussed, (tilde{f}) does not introduce relevant poles, and the forced solution has a form similar to the unforced one.
The harmonic model is, in fact, the minimal model sharing characteristics with our observed structures, and the correspondence between Eqs. (12) and  (13) will allow to estimate the oscillation features. In the limit of small (zeta), (tilde{R}_H(k)) has a sharp peak, characterized by a large quality factor (Q equiv k^star /Delta k), where (Delta k) is the bandwidth at half-height of (tilde{R}(k)) around (k^star)43. First, we see that the position of the peak of (tilde{R}_H) approximately gives the oscillation mode (kappa), according to (k^star = k_0sqrt{1-2zeta ^2} = kappa + mathcal{O}(zeta ^2)). Second, the bandwidth is related to the decay-length through (Delta k = 2/bar{x} + mathcal{O}(zeta ^2))44.
Putting all together, as long as (tilde{R}(k)) resembles the bell-shaped form of (tilde{R}_H(k)), we can use the following estimates, which are correct for the harmonic case to first order in (zeta):

$$begin{aligned} bar{k} simeq underset{k}{text {arg max}},(tilde{R}) equiv k^star quad text { and }quad bar{x} simeq frac{2}{Delta k} ,. end{aligned}$$
(14)

The expression for (bar{x}) is also valid in the overdamped limit (large (zeta) in the harmonic model), in which case the maximum is located at (k^star =0).
The adequacy of the harmonic framework as an approximation to the response function of our model, (tilde{R}(k)), is illustrated in Fig. 6. In the case (D=2times 10^{-1}), the harmonic response is able to emulate (tilde{R}(k)). Then, if the harmonic approximation holds, one expects that the estimates given by Eq. (14) should work for the population dynamics case. In fact, they do work, as we will see below. Differently, when (D=2times 10^{-4}), (tilde{R}(k)) does not follow the harmonic shape, it is multipeaked and the dominant mode observed in the simulations is not given by the absolute maximum.
Figure 6

Comparison of (tilde{R}(k)) with the harmonic response (tilde{R}_H(k)), both normalized to their maximal values. (tilde{R}(k)) of our model, given by Eq. (11) (solid lines) and harmonic response (tilde{R}_H(k)), given by Eq. (12) (dashed lines), where the values of (k_0) and (zeta) were obtained by fitting Eq. (12) to (tilde{R}(k)). In all cases, (q=0.5), (ell =2) and two different values of D shown in the legend were considered. Notice that for (D=2times 10^{-1}), the response can be described by the harmonic approximation. For (D=2times 10^{-4}), the response is multipeaked, indicating that the harmonic approximation fails. In fact, the dominant mode observed in the simulations is not given by the absolute maximum, but by the small hump at (ksimeq 2.1), as predicted by the analysis of complex poles.

Full size image

In Fig. 5, we compare the values of (bar{k}) and (bar{x}) extracted from the numerical solutions of Eq. (3) with those estimated by Eq. (14) (dashed lines) and, more accurately, with those predicted from the poles of (tilde{R}(k)) (solid lines), which perfectly follow the numerical results. The harmonic estimates are shown in the full abscissa ranges, as a reference, even in regions where the approximation is not expected to hold, because discrepancies give an idea of the departure from the harmonic response.
Figure 5c shows outcomes for a fixed (q 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