Ecology

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

1.
Luginbill, P. The fall armyworm. USDA Tech. Bull. 34, 1–91 (1928).
Google Scholar 
2.
Chandrasena, D. I. et al. Characterization of field-evolved resistance to Bacillus thuringiensis-derived Cry1F delta-endotoxin in Spodopterafrugiperda populations from Argentina. Pest Manag. Sci. 74, 746–754. https://doi.org/10.1002/ps.4776 (2018).
CAS  Article  PubMed  Google Scholar 

3.
Farias, J. R. et al. Frequency of Cry1F resistance alleles in Spodoptera frugiperda (Lepidoptera: Noctuidae) in Brazil. Pest Manag. Sci. 72, 2295–2302. https://doi.org/10.1002/ps.4274 (2016).
CAS  Article  PubMed  Google Scholar 

4.
Farias, J. R. et al. Field-evolved resistance to Cry1F maize by Spodoptera frugiperda (Lepidoptera: Noctuidae) in Brazil. Crop Prot. 64, 150–158 (2014).
Article  ADS  Google Scholar 

5.
Huang, F. et al. Cry1F resistance in fall armyworm Spodoptera frugiperda: Single gene versus pyramided Bt maize. PLoS ONE 9, e112958. https://doi.org/10.1371/journal.pone.0112958 (2014).
CAS  Article  PubMed  PubMed Central  ADS  Google Scholar 

6.
Storer, N. P. et al. Discovery and characterization of field resistance to Bt maize: Spodoptera frugiperda (Lepidoptera:Noctuidae) in Puerto Rico. J. Econ. Entomol. 103, 1031–1038. https://doi.org/10.1603/Ec10040 (2010).
Article  PubMed  Google Scholar 

7.
Ganiger, P. C. et al. Occurrence of the new invasive pest, fall armyworm, Spodoptera frugiperda (JE Smith) (Lepidoptera: Noctuidae), in the maize fields of Karnataka, India. Curr. Sci. India 115, 621–623 (2018).
CAS  Article  Google Scholar 

8.
Nagoshi, R. N. Evidence that a major subpopulation of fall armyworm found in the Western Hemisphere is rare or absent in Africa, which may limit the range of crops at risk of infestation. PLoS ONE 14, e0208966. https://doi.org/10.1371/journal.pone.0208966 (2019).
CAS  Article  PubMed  PubMed Central  Google Scholar 

9.
Nagoshi, R. N. et al. Southeastern Asia fall armyworms are closely related to populations in Africa and India, consistent with common origin and recent migration. Sci. Rep. 10, 1421. https://doi.org/10.1038/s41598-020-58249-3 (2020).
CAS  Article  PubMed  PubMed Central  ADS  Google Scholar 

10.
Shylesha, A. N. et al. Studies on new invasive pest Spodopterafrugiperda (J. E. Smith) (Lepidoptera: Noctuidae) and its natural enemies. J. Biol. Control 32, 145–151. https://doi.org/10.18311/jbc/2018/21707 (2018).
Article  Google Scholar 

11.
DPIRD, G. o. W. A. Fall armyworm in Western Australia. http://www.agric.wa.gov.au/plant-biosecurity/fall-armyworm-western-australia (2020).

12.
Pair, S. D. & Sparks, A. N. in Long-range migration of moths of agronomic importance to the United States and Canada: Specific examples of occurrence and synoptic weather patterns conducive to migration (ESA Symposium, 1982). Vol. ARS-43 (ed A. N. Sparks) 25–33 (USDA Miscellaneous Publication, 1986).

13.
Mitchell, E. R. et al. Seasonal periodicity of fall armyworm, (Lepidoptera, Noctuidae) in the Caribbean basin and northward to Canada. J. Entomol. Sci. 26, 39–50 (1991).
Article  Google Scholar 

14.
Nagoshi, R. N., Meagher, R. L. & Hay-Roe, M. Inferring the annual migration patterns of fall armyworm (Lepidoptera: Noctuidae) in the United States from mitochondrial haplotypes. Ecol. Evol. 2, 1458–1467 (2012).
Article  Google Scholar 

15.
Westbrook, J. K. Noctuid migration in Texas within the nocturnal aeroecological boundary layer. Integr. Comp. Biol. 48, 99–106 (2008).
Article  Google Scholar 

16.
16Danthanarayana, W. in Proceedings in life sciences (ed International Congress of Entomology) (Springer, Hamburg, 1986).

17.
Westbrook, J. K., Nagoshi, R. N., Meagher, R. L., Fleischer, S. J. & Jairam, S. Modeling seasonal migration of fall armyworm moths. Int. J. Biometeorol. 60, 255–267 (2016).
CAS  Article  ADS  Google Scholar 

18.
Pashley, D. P. The current status of fall armyworm host strains. Fla Entomol. 71, 227–234 (1988).
Article  Google Scholar 

19.
19Pashley, D. P. in Electrophoretic Studies on Agricultural Pests (eds H. D. Loxdale & J. der Hollander) 103–114 (Oxford University Press, Oxford, 1989).

20.
Juárez, M. L. et al. Host association of Spodoptera frugiperda (Lepidoptera: Noctuidae) corn and rice strains in Argentina, Brazil, and Paraguay. J. Econ. Entomol. 105, 573–582. https://doi.org/10.1603/Ec11184 (2012).
Article  PubMed  Google Scholar 

21.
Murúa, M. G. et al. Demonstration using field collections that Argentina fall armyworm populations exhibit strain-specific host plant preferences. J. Econ. Entomol. 108, 2305–2315 (2015).
Article  Google Scholar 

22.
Nagoshi, R. N., Silvie, P., Meagher, R. L., Lopez, J. & Machados, V. Identification and comparison of fall armyworm (Lepidoptera: Noctuidae) host strains in Brazil, Texas, and Florida. Ann. Entomol. Soc. Am. 100, 394–402 (2007).
CAS  Article  Google Scholar 

23.
Levy, H. C., Garcia-Maruniak, A. & Maruniak, J. E. Strain identification of Spodoptera frugiperda (Lepidoptera: Noctuidae) insects and cell line: PCR-RFLP of cytochrome oxidase C subunit I gene. Fla Entomol. 85, 186–190 (2002).
CAS  Article  Google Scholar 

24.
Nagoshi, R. N. The fall armyworm triose phosphate isomerase (Tpi) gene as a marker of strain identity and interstrain mating. Ann. Entomol. Soc. Am. 103, 283–292. https://doi.org/10.1603/An09046 (2010).
CAS  Article  Google Scholar 

25.
Nagoshi, R. N., Silvie, P. & Meagher, R. L. Comparison of haplotype frequencies differentiate fall armyworm (Lepidoptera: Noctuidae) corn-strain populations from Florida and Brazil. J. Econ. Entomol. 100, 954–961 (2007).
Article  Google Scholar 

26.
Nagoshi, R. N., Meagher, R. L. & Jenkins, D. A. Puerto Rico fall armyworm has only limited interactions with those from Brazil or Texas but could have substantial exchanges with Florida populations. J. Econ. Entomol. 103, 360–367 (2010).
Article  Google Scholar 

27.
Nagoshi, R. N. et al. Genetic characterization of fall armyworm infesting South Africa and India indicate recent introduction from a common source population. PLoS ONE 14, e021775 (2019).
Google Scholar 

28.
Nagoshi, R. N., Fleischer, S. & Meagher, R. L. Demonstration and quantification of restricted mating between fall armyworm host strains in field collections by SNP comparisons. J. Econ. Entomol. 110, 2568–2575 (2017).
CAS  Article  Google Scholar 

29.
Nagoshi, R. N., Goergen, G., Du Plessis, H., van den Berg, J. & Meagher, R. Genetic comparisons of fall armyworm populations from 11 countries spanning sub-Saharan Africa provide insights into strain composition and migratory behaviors. Sci. Rep. UK 9, 8311 (2019).
Article  ADS  Google Scholar 

30.
Nagoshi, R. N. & Meagher, R. L. Using intron sequence comparisons in the triose-phosphate isomerase gene to study the divergence of the fall armyworm host strains. Insect Mol. Biol. 25, 324–337. https://doi.org/10.1111/imb.12223 (2016).
CAS  Article  PubMed  Google Scholar 

31.
Meagher, R. L. & Nagoshi, R. N. Population dynamics and occurrence of Spodoptera frugiperda host strains in southern Florida. Ecol. Entomol. 29, 614–620 (2004).
Article  Google Scholar 

32.
Nagoshi, R. N. et al. Genetic characterization of fall armyworm (Lepidoptera: Noctuidae) host strains in Argentina. J. Econ. Entomol. 105, 418–428. https://doi.org/10.1603/Ec11332 (2012).
CAS  Article  PubMed  Google Scholar 

33.
Meagher, R. L. & Nagoshi, R. N. Differential feeding of fall armyworm (Lepidoptera: Noctuidae) host strains on meridic and natural diets. Ann. Entomol. Soc. Am. 105, 462–470. https://doi.org/10.1603/An11158 (2012).
Article  Google Scholar 

34.
Meagher, R. L., Nagoshi, R. N., Stuhl, C. & Mitchell, E. R. Larval development of fall armyworm (Lepidoptera: Noctuidae) on different cover crop plants. Fla Entomol. 87, 454–460 (2004).
Article  Google Scholar 

35.
Prowell, D. P., McMichael, M. & Silvain, J. F. Multilocus genetic analysis of host use, introgression, and speciation in host strains of fall armyworm (Lepidoptera: Noctuidae). Ann. Entomol. Soc. Am. 97, 1034–1044 (2004).
CAS  Article  Google Scholar 

36.
Nagoshi, R. N. et al. Genetic characterization of fall armyworm (Spodoptera frugiperda) in Ecuador and comparisons with regional populations identify likely migratory relationships. PLoS ONE 14, e0222332. https://doi.org/10.1371/journal.pone.0222332 (2019).
CAS  Article  PubMed  PubMed Central  Google Scholar 

37.
Nagoshi, R. N. et al. Haplotype profile comparisons between Spodoptera frugiperda (Lepidoptera: Noctuidae) populations from Mexico with those from Puerto Rico, South America, and the United States and their implications to migratory behavior. J. Econ. Entomol. 108, 135–144 (2015).
CAS  Article  Google Scholar 

38.
Nagoshi, R. N. et al. Fall armyworm migration across the Lesser Antilles and the potential for genetic exchanges between North and South American populations. PLoS ONE 12, e0171743 (2017).
CAS  Article  Google Scholar 

39.
Peck, S. B., Heraty, J., Landry, B. & Sinclair, B. J. Introduced insect fauna of an oceanic archipelago: The Galápagos Islands, Ecuador. Am. Entomol. 44, 218–237. https://doi.org/10.1093/ae/44.4.218 (1998).
Article  Google Scholar 

40.
Zapata, F. & Granja, M. M. Optimizing marine transport of food products to Galapagos: advances in the implementation plan. http://www.galapagos.org/wp-content/uploads/2012/04/trans1-optimizing-marine-transport.pdf (2009–2010).

41.
Toral-Granda, M. V. et al. Alien species pathways to the Galapagos Islands, Ecuador. PLoS ONE 12, e0184379. https://doi.org/10.1371/journal.pone.0184379 (2017).
CAS  Article  PubMed  PubMed Central  Google Scholar 

42.
Nagoshi, R. N. et al. The genetic characterization of fall armyworm populations in Ecuador and its implications to migration and pest management in the northern regions of South America. PLoS ONE 15, e0236759. https://doi.org/10.1371/journal.pone.0236759 (2020).
CAS  Article  PubMed  PubMed Central  Google Scholar 

43.
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).
Article  Google Scholar 

44.
Saitou, N. & Nei, M. The neighbor-joining method—A new method for reconstructing phylogenetic trees. Mol. Biol. Evol. 4, 406–425 (1987).
CAS  PubMed  Google Scholar 

45.
Clements, M. J., Kleinschmidt, C. E., Maragos, C. M., Pataky, J. K. & White, D. G. Evaluation of inoculation techniques for fusarium ear rot and fumonisin contamination of corn. Plant Dis. 87, 147–153 (2003).
CAS  Article  Google Scholar 

46.
Leigh, J. W. & Bryant, D. POPART: Full-feature software for haplotype network construction. Methods Ecol. Evol. 6, 1110–1116 (2015).
Article  Google Scholar 

47.
Murúa, G. M. et al. Fitness and mating compatibility of Spodoptera frugiperda (Lepidoptera: Noctuidae) populations from different host plant species and regions in Argentina. Ann. Entomol. Soc. Am. 101, 639–649 (2008).
Article  Google Scholar 

48.
Librado, P. & Rozas, J. DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452 (2009).
CAS  Article  Google Scholar 

49.
Stein, A. F. et al. NOAA’s HYSPLIT atmospheric transport and dispersion modeling system. Bull. Am. Meteorol. Soc. 96, 2059–2077. https://doi.org/10.1175/Bams-D-14-00110.1 (2015).
Article  ADS  Google Scholar  More

Physiological responses
Small fragments from five source colonies from the two experimental sites (N- and O-sites) were used to conduct a reciprocal transplant experiment in Maunalua Bay, Hawaii (Fig. 1). The results revealed clear physiological response differences between the two populations. The transplantation resulted in a significant reduction in the average tissue layer thickness (TLT) in only one treatment: O-corals transplanted to N-site (O → N) (Tukey-HSD, P-adj  2 at FDR = 0.01. Proteins associated with key GO terms were colored in different colors, and the top 10 abundant proteins in each population are annotated. The bottom bars indicate the total numbers of significantly abundant proteins for each population.

Full size image

Response difference in transplant to the offshore site (N → O vs. O → O)
A total of 3236 distinct coral proteins were identified at O-site: 2217 (68.5%) were shared between the two populations, 656 unique to N → O corals, and 363 to O → O corals (Fig. S1C). GO analysis identified 35 enriched terms specific to N → O, which involved amino acid biosynthetic process, ATP metabolic process, TCA cycles, fatty acid oxidation, and monosaccharide metabolic process. There were 15 specific GO terms in O → O corals, including nucleotide monophosphate biosynthetic process, intracellular protein transport, vesicle organization, and GTP binding (SI.2B).
Quantitative analysis on protein abundances indicated a total of 665 proteins to be significantly differentially abundant at O-site: N → O corals had 155 abundant-proteins, and O → O corals had 510 abundant-proteins (Fig. 3B). GO analysis resulted in identifying 39 enriched terms from abundant proteins in O → O corals, while only one met the cutoff in N → O corals (SI.2B). Although the number of abundant-proteins and enriched terms identified in O → O corals were relatively high, the enriched terms predominantly consisted of cellular functions related to protein translation; organonitrogen biosynthetic process and organic acid metabolic process, both leading to single child terms for BP, CC, and MF (tRNA aminoacylation for protein translation, cytosolic large ribosomal subunit, and tRNA aminoacyl ligase activity). The enriched term in N → O corals was a non-specific term of ‘extracellular region’, indicating that despite the higher number of abundant-proteins, the main functional difference between N → O and O → O corals was an enhanced protein translation activity in O → O corals.
Response comparisons to cross transplantation
Effects of cross transplantation yielded a more diverse proteomic stress-response in O-corals as they moved nearshore than N-corals as they were moved offshore (Fig. S2). The total number of abundant-proteins between the sites was much higher for O-corals (440, O → N vs. O → O) than N-corals (135, N → N vs. N → O) (Table S1), and the number of unique GO terms identified between the sites was also higher in O-corals (69, SI.2C) than in N-corals (46, SI.2D). The number of overlapping proteins between the sites was lower in O-corals than in N-corals (70% vs. 79%), and log-fold changes of all identified proteins between the sites were significantly larger for O-corals than N-corals (Wilcoxon Rank-Sum test, P = 6.02 × 10–9), all emphasizing the larger metabolic reshuffling needed to respond to cross transplantation in O-corals. GO enrichment analysis indicated that N-corals responded to transplantation to O-site with increased abundance of proteins involved in amino acid biosynthesis, fatty acid beta oxidation, TCA cycle, chitin catabolism, coenzyme biosynthesis and translational initiation. O-corals responded to transplantation to N-site by increasing the abundance of proteins associated with detoxification, antioxidant activity, protein complex subunit organization, and multiple metabolic processes (amino acid, fatty acid, ATP, monosaccharide, and carbohydrate derivative) (SI.2E). The shared responses between the cross-transplanted corals (N → O and O → N corals) included increased proteins involved in fatty-acid beta oxidation, TCA cycle, carbohydrate derivative catabolic process, pyridoxal phosphate binding, and ‘oxidoreductase activity acting on the CH-CH group of donors with flavin as acceptor’, likely representing the effects of transplantation to a non-native environment.
Proteome patterns across the four treatments
Comparing enriched GO terms across all treatments (SI.2E) highlighted the unique state of O → N corals; O → N corals had a much higher number of uniquely enriched GO terms (n = 27) compared to those in the rests (4 in O → O, 5 in N → N, and 15 in N → O corals). The most notable difference among the treatments was enrichment of detoxification and antioxidant activity exclusively in O → N corals (Fig. 4). Also, lipid oxidation was highly enriched in O → N corals with four terms associated to this category identified (Fig. 4, SI.2E).
Figure 4

Enriched GO terms uniquely identified to specific treatment groups. Treatment groups are shown in the right column (e.g. N-coral = N-corals at both sites, N-site = N- and O-corals at N-site, CrossT = cross transplantation). The heat-map represents P-values for the associated GO terms. The GO terms are grouped by the parent–child terms with the most parent term in bold (for values, see SI-2E).

Full size image

Examining the relative abundance of individual proteins associated with detoxification (‘detox-proteins’) revealed the following interesting patterns. (1) Distinct sets of proteins were abundant in different treatments, rather than all detox-proteins to be elevated in one treatment, and the direction and magnitude of responses to transplantation were protein specific and varied between populations (Fig. S4A). (2) Two peroxiredoxin (Prx) proteins, Prx-1 (m.6147) and Prx-6 (m.9595), dominated the relative abundance of detox-proteins by having over an order of magnitude higher abundance values, and they were consistently more abundant in N-corals than O-corals (ave. 44%, Kruskal Test, P = 0.004–0.01) (Fig. S4B, SI.1B). (3) Some proteins with the same or similar annotations had contrasting responses between the populations. For example, Prx-4 (m.17739), which belongs to the same subfamily as Prx-1, was significantly more abundant in O-corals at both sites (Fig. S4B, SI.2F,G), while Prx-1 was more abundant in N-corals. Similarly, seven peroxidasin (PXDN) homologs were identified, of which m.17686 was significantly more abundant in O → N corals, while m.9432 was significantly more abundant in N → N corals (Fig. S4B, SI.2F), suggesting that the two populations potentially utilize different class/kind of enzymes as primary proteins in detoxification/antioxidant pathways. Of the seven PXDN homologs, two (m.1440, m.9432) were consistently higher in N-corals, two (m.10928, m.15200) were consistently higher in O-corals, and three (m.12572, m.17686, m.9657) increased abundance at N-site in both corals, but m.12572 and m.17686 being higher in O-corals, while m.9657 higher in N-corals (Fig. S3B).
To ascertain that the proteins with the same annotations are indeed different proteins, sequences of matched peptides were assessed for those that showed contrasting responses. The pairwise comparison of Prx-1 and Prx-4 showed only seven of the total 65 peptides (11%) were identical between the two, revealing that these protein sequences are significantly different and they each have unique peptides that be detected and quantified accurately (SI.1C1). Similarly the majority of PXDN-like proteins identified had no overlapping peptides between the contrasting pairs (0–19%, median = 0, SI.1C2), indicating that corals possess multiple types of PXDN, and N- and O-corals respond to stressors with different sets of PXDN.
In addition to lipid oxidation being significantly enriched in O → N corals, a single term (fatty acid beta-oxidation,) was also enriched in N → O corals, which suggests that cross-transplantation had an effect on lipid oxidation processes. However, the abundances of most proteins associated with lipid oxidation were higher in O-corals than N-corals at both sites (Fig. S4A). Statistically, three proteins (medium-chain sp acyl-CoA:m.22274, very-long-chain sp. acyl-CoA:m.17984, and trifunctional enzyme subunit alpha:m.6724) showed a difference in abundance between the two populations at N-site (Fig. S4C) and one (isovaleryl-CoA dehydrogenase:m.27714) at O-site, all of which were higher in O-corals than N-corals. More