Philippa Kaur

SubjectsThe study comprised 34 men (19–26 years, mean 23 years; body mass index 19–31 kg/m2, mean 24 kg/m2) with no physical disabilities or mental disorders, including color blindness and claustrophobia30,31. All subjects were informed of the aims, the study procedure, and the financial compensation for participation, and were asked to follow the rules of the study. Prior to each experiment, subjects underwent short-term starvation31,54 (18–20 h; no food except pure water after lunch (12–1 pm) or dinner (6–7 pm), no later than 1 pm or 7 pm, respectively, one the day before the test), no medical treatments, and normal sleep (at least 6 h, between 10 pm the day before the test day to 8 am on the test day)31. Prior to starting each experiment, subjects were stabilized on a chair for ~ 10 min in a room next to the testing room. Based on an assessment with a pre-experiment questionnaire and the first blood glucose level, measured before starting the experiment (see “Geomagnetic orientation assay” section below), subjects who had not followed these rules were not allowed to take the test on the day and the test was postponed. The study was approved by the Institutional Review Board of Kyungpook National University and all the procedures followed the regulations for human subject research. Informed consent was obtained from all subjects.Modulation of GMFThe ambient GMF in the testing room had a total intensity 45.0 μT, inclination 53°, and declination − 7° (Daegu city, Republic of Korea); the total intensity of 50.0 μT in our previous study31 was changed due to a reconstruction of the building; 45.0 μT was maintained throughout the period of this study. To provide the subjects with various GMF-like magnetic fields (i.e., by modulating of total intensity, inclination, or direction of magnetic north), the coil system from our previous studies6,7,31 was used. It comprised three double-wrapped, orthogonal, rectangular Helmholtz coils (1.890 × 1.890 m, 1.890 × 1.800 m, and 1.980 × 1.980 m for the north–south, east–west, and vertical axes, respectively), electrically-grounded with copper mesh shielding. The subject was seated on a rotatable plastic chair with no metal components, at the center of the three-dimensional coils with his head positioned in the middle space of the vertical axis of the coils. To modulate the geomagnetic north, each pair of coils was supplied with direct current from a power supply (MK3003P; MKPOWER, Republic of Korea). The magnetic field was measured using a 3-axis magnetometer (MGM 3AXIS; ALPHALAB, USA); the field homogeneity at the position of the subject’s head was found to be 95%. The testing room was shielded by a six-sided Faraday cage comprising 10 mm thick aluminum plates, and was grounded during the entire experiment40. Background electromagnetic noise was measured inside the coils at the start and the end of each experimental day. It was attenuated by the Faraday cage more than 200-fold over the range from 500 Hz to 100 MHz as described in detail in our previous study31. The 60 Hz power frequency magnetic field was no more than 2 nT (3D NF Analyzer NFA 1000; Gigahertz Solutions, Germany). All electronic devices were placed outside the Faraday cage during the experiments, with the exception of the switch button module for GMF modulation and the antenna for generating the oscillating magnetic fields. The temperature experienced by the subjects was maintained at 25 ± 0.5 °C (Data logger 98,581; MIC Meter Industrial, Taiwan) by air circulation through the honeycomb on the ceiling of the Faraday cage31.Geomagnetic orientation assayAdopting a two-alternative forced choice (2-AFC) paradigm33,34, a geomagnetic orientation assay was conducted similar to our previous study31. Experiments were performed at 09:30–11:30 am or 1:00–5:00 pm (local time, UTC + 09:00) (each experiment: 50 min–1 h 10 min; mean ≈ 1 h, which was shorter by approximately 30 min than that in the previous study: 1 h 20 min–1 h 40 min; mean ≈ 1 h 30 min). Depending on the experiment, starved or unstarved subjects were tested individually. Prior to each experiment, the subjects were asked to remain with their heads facing the front, with eyes closed and earmuffs on during the experiment. In particular, they were asked to concentrate on sensing, if they could, the ambient geomagnetic north during the association phase, and to use the sensed information, depending on the experiment, to orient toward one of the two modulated magnetic norths (0°/180° for magnetic north–south axis or 90°/270° for magnetic east–west axis, rotated clockwise with respect to the ambient geomagnetic north) during the test phase. Subjects were instructed to avoid distracting thoughts and to think immediately “which direction is modulated magnetic north?” whenever they were distracted during the test phase, or felt they were being biased by experiences from earlier experiments. While seated on the rotatable chair, the subject’s blood glucose level was measured shortly before the first session and immediately after each session with eyes open except in the ‘dark’ experiment (Accu-Chek Guide; Roche, Germany)31. If the determined value before the first session varied by more than 15% relative to the mean (Table S2)31, the experiment was postponed and repeated at a later date (approximately 2% of experiments). The subjects were stabilized with eyes closed for 2 min before the first trial in the absence of visual, auditory, olfactory, and haptic sensory cues. For the ‘dark’ experiment (light intensity ≈ 0 lx), subjects wore home-made ‘blind’ goggles and were stabilized with eyes closed for 5 min55,56, and then asked whether they could see any light. If they could, the goggles were adjusted to prevent leakage of light, and the subject then had another 5 min of stabilization with eyes closed before starting the experiment. The subjects were illuminated with light from a filtered/non-filtered diffused light-emitting diode, depending on the experiment (Table S1). The home-made filter goggles were worn throughout the experiment, including the association and test phase, when required. The goggles contained glass filters (Tae Young Optics, Republic of Korea) to provide the eyes with particular wavelengths of light (Spectrometer USB4000-UV-VIS, Ocean Optics, USA) (Fig. S1). Each experiment consisted of 16 sequential trials for ‘no-association’ and ‘food-association’. For the food-association, a subject facing toward the ambient geomagnetic north was gently provided with a chocolate chip31 on his right palm by an experimenter, and given 30 s to eat it, while during no-association trials, food was not provided during the association phase. After a subsequent 5 s interval, the experimenter gently touched the subject’s right thenar area using a paper rod, as the cue to start the test. One of the two modulated magnetic north directions, depending on the experiment, was randomly provided 3 s before the cue for the test. Each of the modulated magnetic north directions was provided eight times for the no-association and food-association sessions. Subjects were informed of the nearly equal probability for each of the modulated magnetic north directions before each experiment. With the touch cue, subjects were asked to rotate freely toward any direction (clockwise or counterclockwise) by themselves (1–4 cycles of two-thirds rotation) and try to sense the direction of the modulated magnetic north during a 1 min period. Rotation was allowed within the rotation angle (− 30° to 210° for the magnetic north–south axis or − 120° to 120° for the east–west axis, depending on experiments, with respect to the ambient magnetic north), which was confined by the plastic stool (Fig. 1A) touching the left or right ankle of the subjects. When subjects determined the direction of the magnetic north, they stopped rotating to face toward the direction and lifted their right hand to indicate the direction to the experimenter. The direction was measured by the experimenter at 10° intervals using the scale on the walls of the Faraday cage31. A prerequisite for correct orientation was that the subject indicated the direction within the range of 30° to the both sides with respect to the magnetic cardinal directions, which was instructed to the subjects before each experiment. When the direction the subjects indicated was out of the 30° range, the trial was not included in the data and was repeated (approximately 0.63% of trials). Before the subsequent trial, the subject was gently rotated to face toward the ambient geomagnetic north and then rested for 5 s. For the ‘dark’ experiment, subjects were asked whether they could see any leaked light immediately after the last measurement of blood glucose level at the end of experiment. If the subject could see leaked light, the experiment was nullified and repeated later on (approximately 3% of experiments; 2/68). All experiments were performed in a double-blind fashion. The experimenter who conducted the orientation assay knew whether a subject was starved or not, wearing filter goggles, and food-associated or not, but did not know the random magnetic north sequences that were controlled by the personal computer (PC) system. Another experimenter responsible for analyzing the data did not know whether the subject was starved or not, the experimental conditions, including light wavelengths, or whether an oscillating magnetic field had been provided to the subjects. Thus, none of the experimenters were aware of all the subject information and data during the experiments and data analysis. The correct orientation rate was calculated by (the number of correct orientation trials/total number of trials) (raw data, Appendix S3). All the subjects participated in all the experiments performed in random order with an interval of at least 3 days between experiments. After each experiment, the subjects were asked to answer a post-experiment questionnaire about whether they closed their eyes when required during the entire period of the experiment. In cases when a subject did not maintain closed eyes, the experiment was repeated (approximately 1% of experiments). For each subject, a preliminary experiment on the “magnetic north–south axis” was conducted twice (unstarved and starved for each) with no goggles for adaption to the experimental procedure. These data were not included in the results.Experiments with oscillating magnetic fieldsExperiments with oscillating magnetic fields were performed using the standard geomagnetic orientation assay described above. To produce the oscillating magnetic fields, oscillating currents from a function generator (AFG3021; Tektronix, USA. For each magnetic field, sweep of 500 ms; interval of 1 ms. See Fig. S6A) were amplified (ENI 2100L RF power amplifier; Bell Electronics, USA) and fed into a calibrated coil antenna (30 cm diameter, 6509 loop antenna; ETS-LINDGREN, USA) mounted on a wooden frame, comprised of a single winding of coaxial cable. The oscillating magnetic fields were measured daily, before the first and after the last experiment of the day, using a spectrum analyzer (SPA-921TG; Com-Power, USA) with a calibrated loop antenna (48 cm diameter, AL-130R; Com-Power, USA) and a calibrated magnetometer (Probe HF 3061, NBM-550; Narda, Germany). Magnetic field intensities were measured on the glabella of the subjects; variations in intensity between subjects due to different seating heights were less than 10% of the average values (Table S3). The function generator and amplifier were placed outside the Faraday cage, and switched on during the dummy load control experiments with no signal from the PC system. The band widths of the monochromatic magnetic fields, i.e., 1.260 and 1.890 MHz were 0.020 and 0.019 MHz (“average”, √3 kHz), respectively, at the bottoms of the peaks. During the test phase, the maximum values of magnetic noise on the glabella of subjects including the dummy load did not exceed the following values: (1) 5 Hz–9 kHz; 2 nT/√ 2 kHz of “average” and 8 nT/√ 9 kHz of “max-hold” (0.05 nT/√ 2 kHz of “average” and 5 nT/√ 9 kHz of “max-hold” in the dummy load) (3D NF Analyzer NFA 1000; Gigahertz Solutions, Germany); (2) 9 kHz–500 kHz; 5 nT/√ 3 kHz of “average” and 8 nT/√ 3 kHz of “max-hold” (≈ 0 nT/√ 3 kHz of “average” and ≈ 1 nT/√ 3 kHz of “max-hold” in the dummy load) (the AL-130R antenna) (Fig. S6C); and (3) 500 kHz–30 MHz; 0.006 nT of 3.780 MHz harmonic in the 1.260 MHz, 0.03 nT of 5.670 MHz harmonic in the 1.890 MHz, and ≈ 0 nT in the dummy load (/√ 10 kHz of “average”) (Fig. S6B), and 0.15 nT/√ 10 kHz of “max-hold” at the same frequencies above and ≈ 0 nT in the dummy load (the AL-130R antenna).Statistical analysisTo determine the significance of orientation data from the 2-AFC paradigm, a one-sample t-test (test mean: 0.5), paired sample t-test, or two-sample t-test was performed using Origin software 11 (Origin, USA). To verify the reasonability of the t-tests, all data sets were checked using the Anderson–Darling test if the data follow a normal distribution (Appendix S4). To evaluate the difference between the means of two data sets when at least one of them did not show a normal distribution, the percentile bootstrap method57 was used (95% confidence interval, see Fig. S2, Appendices S1 and S2 for raw data). To analyze the blood glucose data, a paired sample t-test was used. Based on the results of previous study31, to describe different response groups of magnetic orientation in the 2-AFC paradigm, a principal component analysis36,37 was conducted on correct orientation rates by starved subjects, with no association/food-association under the full wavelength or  > 400 nm light conditions using SPSS 23 (IBM, USA). Following the principal component analysis calculation, the k-means clustering algorithm—one of the unsupervised learning methods—was used to objectively classify the groups58. The number of groups was two, and the distance between the center of the cluster and all points was Euclidean distance. The classification boundary was marked with the perpendicular bisector from the centers of the two groups. The first two principal components accounted for a significant portion of the total variance (73.1%; PC1 = 40.8%, PC2 = 32.3%). Statistical values are presented as mean ± SEM. More

Stability analysis of the deterministic modelSolving (left( x(t) times left( r_{p}x(t)^{(alpha )}left( 1-big (frac{x(t)}{K}big )^{beta }right) – nx(t)^{(delta )} right) right) =0), we obtain two stable and one unstable equilibrium points for the model. One stable equilibrium is trivial, i.e., (x(t)=0), another stable equilibrium point being the non-zero satisfying (left( r_{p}x(t)^{(alpha )}left( 1-big (frac{x(t)}{K}big )^{beta }right) – nx(t)^{(delta )} right) =0). Figure 1a shows three different equilibrium points of the model. In addition to the equilibrium, the model has two inflection points (Fig. 1a). At these inflection points the absolute growth rates are minimum and maximum. The density vs relative proliferation rate (RPR) profile of the model shows that the model can attain negative RPR for a positive cell density, suggesting that the model can portray the Allee phenomenon (Fig. 1b). Figure 1c,d portray the proliferation and decay phases, respectively through the model.Figure 1Growth dynamics of the proposed model: (a) Absolute proliferation rate (APR) profile considering (r_{p}=0.13), (K=1.43), (n=0.0095), (alpha =1.15), (beta =0.99) and (delta =0.2); (b) RPR profiles for different n and other same constant model parameters; (c) Cell population survive for (r_{p}=0.13), (K=1.43), (n=0.0095), (alpha =1.15), (beta =0.99) and (delta =0.2) with the initial cell density 0.1; (d) The population goes to extinction for the initial cell density 0.06 with the same constant parameters.Full size imageThe solution of the deterministic model finally provides two theorems.
Theorem 1

(x^{*}approx K -Kleft( frac{Big (beta r_{p}K^{alpha }+n delta K^{delta }Big )-sqrt{Big (beta r_{p}K^{alpha }+n delta K^{delta }Big )^{2}-2 left( 2 alpha beta r_{p}K^{alpha } +beta (beta -1)r_{p}K^{alpha }+delta (delta -1)nK^{delta } right) nK^{delta }}}{left( 2 alpha beta r_{p}K^{alpha } +beta (beta -1)r_{p}K^{alpha }+delta (delta -1)nK^{delta } right) }right)) is the conditional MSSCD for the intercellular-interaction-induced proliferative cells. The conditional threshold density for cell-proliferation upon interaction is (x^{*}=K -Kleft( frac{Big (beta r_{p}K^{alpha }+n delta K^{delta }Big )+sqrt{Big (beta r_{p}K^{alpha }+n delta K^{delta }Big )^{2}-2 left( 2 alpha beta r_{p}K^{alpha } +beta (beta -1)r_{p}K^{alpha }+delta (delta -1)nK^{delta } right) nK^{delta }}}{left( 2 alpha beta r_{p}K^{alpha } +beta (beta -1)r_{p}K^{alpha }+delta (delta -1)nK^{delta } right) }right)) (proof is in the supplementary information).
Allee and cooperation models are the only extended logistic law other than our model to provide a threshold population size for growth or proliferation. Our proposed model is superior to the Allee and cooperation model as it can detect the conditional threshold cell density for proliferation and regulate the density by its different parameters. For example, One may reduce the conditional threshold density by either regulating the interaction between growth-inhibiting molecules and cells ((delta)) or reducing the inhibiting molecule concentration (n).The conditional MSSCD from Theorem 1 is lower than the carrying capacity of the conventional logistic model due to growth-inhibiting molecules; it provides the expected cell density during culture in a given environment. Theorem 1 also states the set of parameters to control the cell proliferation and get the desired density during such cell cultures. A further question arises knowing this set of parameters: which one of the parameters in the expression is crucial in terms of application purpose? Since the (r_{p}) is the constant proliferation rate for a given cell line, controlling the conditional MSSCD is not possible through (r_{p}). We simulate the distribution of conditional MSSCD for other parametric planes to answer this question. For this, we use the parameter values obtained from the data.

Theorem 2

The RPR is maximum at the cell density (x^{*}= K-Kleft( frac{r_{p}beta K^{alpha -1}+ndelta K^{delta -1}}{2r_{p}alpha beta K^{alpha -1}+r_{p}beta (beta -1)K^{alpha -1}+ndelta (delta -1)K^{delta -1}}right)) for the concave downward profile under the condition (r_{p}alpha (alpha -1){x^{*}}{}^{(alpha -2)}-frac{r_{p}}{K^{beta }}(alpha +beta )(alpha +beta -1){x^{*}}{}^{(alpha +beta -2)}-ndelta (delta -1){x^{*}}{}^{(delta -2)}n) (see the supplementary information). The cell population sustain with any positive initial cell density x(t) and try to stabilize at (x(t)= K(1-frac{n}{r_{p}})^frac{1}{beta }). Therefore, bimodality vanishes and unimodality is observed for the case (alpha =delta) (r_{p} >n). The RPR profile will be concave downward always with the maximum RPR value is at the inflection point (x(t)= K(frac{(r_{p}-n)alpha }{r_{p}(alpha +beta )})^frac{1}{beta }). The deterministic potential function in this case is (U(x)=-Big [(r_{p}-n)frac{x^{(alpha +2)}}{(alpha +2)}-frac{r_{p}}{K^{beta }}frac{x^{(alpha +beta +2)}}{(alpha +beta +2)} Big ]). The minima of this effective potential function will be at (x(t)= K(1-frac{n}{r_{p}})^frac{1}{beta }) which is the maximum stable cell density for (r_{p} >n).
Parameter estimationThe density-RPR and time-density fitting to the scratch assay datasets show a lower RSS for our model than the logistic one for each of the three seeding conditions. The estimated parameters from the RPR fitting through the grid-search are in Table 2. Although the RSS for the RPR fitting of the seeding 2 is very low, the data itself is too scattered in both the upper and lower range for the small cell density. Therefore, there is a chance that regardless of the low RSS value, the fitting for seeding 2 may not reflect the actual estimates of the parameters with the bias in the data set (Fig. 2b). Nevertheless, the density-RPR fittings to the other two seeding density datasets do not suffer from bias (Fig. 2a,c).Table 2 Estimated model parameters from density-RPR fitting of our model.Full size table
Figure 2Our proposed model best fitted the cell density-RPR datasets for all of the seeding conditions generated through the grid-search method.Full size image
Jin et al.1 suggested that their two phase logistic model may share similarities with the Allee effect. However, they did not fit the Allee model stating seeding 2 and 3 were large enough seeding densities. We calculated the conditional threshold density, conditional MSSCD, density at the minimum and maximum RPR for the model from our estimated parameters (Table 3). The conditional threshold cell density calculated from our estimated parameters confirms that the smallest initial seeding density of the dataset was greater than the conditional threshold cell density.Table 3 Calculated cell densities from estimated parameters from our model fitting.Full size tableFigure 3 compares the portrayal of the data through our model with the fitting by Jin et al.1. The blue dashed line is the time-series fitting of the proposed model, and the red-colored line is the time-series fitting of the logistic model to the scratch assay data sets in the Fig. 3. The carrying capacity values are unexpectedly very high in the logistic fit, keeping the model near the exponential phase for the entire dataset. Thus the overall and two phase logistic fits are unrealistic compared to the highest cell density observed in the assay. Also, logistic fitting of the RPR profiles to the data after 18 h does not capture the whole scenario. The green solid and the violet dashed line represent the logistic time-density fit after and before 18 h density profiles respectively. The orange-colored lines in the Fig. 3 are the expected population density as per estimated parameters from the RPR fitting after 18 h data sets. Table 4 enlists all parameters for a comparison between logistic and our model fitting.Figure 3Time series solution of the proposed model and logistic law with comparative RSS for all three seeding conditions.Full size imageTable 4 Logistic model fitting with the Jin et al.1 estimates used in Fig. 3 with the specific colors.Full size tableTrends in cell densities under deterministic set upThe (r_{p}) is fixed for a cell line among all the determining parameters of the conditional MSSCD. n and K vary together with the culture media, flask, and environmental setup. On the other hand, the (alpha), (beta), and (delta) vary together with intercellular-interactions and cellular-interaction with growth-inhibitory molecules, which depend on the medium’s initial cell density per well and fluidity. We observe that the distribution of the conditional MSSCD depends more on the K than the n. There is a chance of overproliferation in the deterministic setup under low n but high K. The cells may die under high n. The cell density at maximum RPR also depends more on K than n (Fig. 4). So the cells should be cultured in the larger flask to achieve maximum proliferativeness.Figure 4The distribution of conditional MSSCD and cell density at maximum RPR in n-K parametric plane.Full size imageThe conditional MSSCD depends more on (beta) than (alpha) (Fig. 5a). The cells may tend to overproliferate under both high (alpha) and (beta). The conditional MSSCD does not exist for a high (delta) and low (beta) depending more on (delta) than (beta). The cells may overproliferate only under a high (beta) and low (delta) (Fig. 5b). The conditional MSSCD also depends more on (delta) than (alpha) showing mostly underproliferation of cells in the (delta ~-alpha) parametric plane. Therefore, the proliferation can be controlled via regulating the interaction between the growth-inhibitory molecules and cells followed by density-regulation through contact-inhibition and cell-cell cooperation (Fig. 5c).Figure 5The distribution of the conditional MSSCD in parametric plane of regulators in the growth law: (a) dependence of the conditional MSSCD on (alpha) and (beta) parameters; (b) dependence of the conditional MSSCD on (delta) and (beta) parameters; (c) dependence of the conditional MSSCD on (alpha) and (delta) parameters.Full size imageThe new cell fitness measure, i.e. cell density at maximum RPR depends more on the (alpha) than the (beta) (Fig. 6a). The cells achieve maximum RPR at a great cell density under the high value of these two parameters. Figure 6b,c suggest that cell density depends only a little on the (delta) under high (alpha) and (beta). Under the low value of these two regulators, a high (delta) always reduces the cell density attaining the maximum RPR, resulting a poor cell-fitness.Figure 6The distribution of cell density at maximum RPR in parametric plane of regulators in the growth law: (a) dependence on (alpha) and (beta) parameters; (b) dependence on (alpha) and (delta) parameters; (c) dependence on (delta) and (beta) parameters.Full size imageStochastic model analysisOur proposed stochastic model (3) can be compared with the general stratonovich stochastic differential equation (frac{dx}{dt}=f(x)+g_{1}(x)epsilon (t)+g_{2}(x)Gamma (t)). Comparing it with our proposed stochastic model we obtain (g_{1}(x)=-x^{delta +1}) and (g_{2}(x)=1). Using the help of47, we get noise induced drift (A(x)=r_{p}x^{alpha +1}left( 1-Big (frac{x}{K}Big )^{beta } right) -nx^{(delta +1)}+D(delta +1)x^{(2delta +1)}-lambda sqrt{DQ}(delta +1)x^{delta }) and noise induced diffusion coefficient (B(x)=Dx^{(2delta +2)}-2lambda sqrt{DQ}x^{(delta +1)}+Q). The cell density at long run can be obtained from the steady state probability density function (SSPDF). The analytical expression of the SSPDF is obtained from the Fokker-Planck equation. The Fokker-Planck equation is (frac{partial P(x, t)}{partial t} =- frac{partial big [ A(x) P(x, t)big ]}{partial x}+ frac{partial ^{2} big [B(x) P(x, t)big ]}{partial x^{2}}), where P(x,t) is the probability density function of the cell population at the time point t. Solving the Fokker-Planck equation we get the SSPDF as (P_{st} (x)= frac{N^{prime }}{B(x)} exp left( int _{x} frac{A(x^{prime })}{B(x^{prime })} dx^{prime }right)) with the normalizing constant (N^{prime }). The value of (N^{prime }) can be obtained from (int _{0}^{infty } P_{st} (x)dx=1).This SSPDF (P_{st} (x)) helps to understand the validity of the proposed stochastic model. Since the number of the data points is too low to fit the stochastic model to the data directly, validation of the stochastic model is challenging in this case. The dataset we used is a time series with 15 data points with three replicates only. An experiment must have many replicates to have a sample with a large sample size so that the SSPDF of cell densities obtained from theoretical findings can be validated with the real observation of cell densities at the steady state. Such datasets with many replicates are rare.So, we generate 2000 sample paths with the help of numerical simulation based on stochastic model 3. We use the parameter values estimated from the fittings of the deterministic model to the seeding condition 1, and we consider some particular values for the two noise intensities and correlation strength ((lambda)) to get a simulated dataset. To achieve the stationary state, we consider sufficiently large time points, and the cell densities at the final time point are used as the data set for the stationary state. We compare the frequency density of cell densities at steady-state of a simulated dataset of 2000 sample paths with the SSPDF obtained from the analytical solution. This comparison shows that the cell density distribution at the steady state matches the steady state probability density function obtained analytically (Fig. 7).In addition, we illustrated the time series generated with the help of stochastic model 3 through numerical technique (Fig. 8). We have plotted the time series data thus obtained for each of the three seeding conditions and in the same figure we also plotted the observed cell densities. The red dots (o) represent the original/experimental dataset of Jin et al.1. The blue dots ((*)) represent the simulated dataset obtained from the stochastic model. This Fig. 8 clarifies our claim that the proposed stochastic model is in good agreement with the actual observation.Figure 7The histogram shows the distribution of cell densities at steady state under additive and multiplicative noises. The blue curve is the SSPDF. The function SSPDF and the distribution of cell densities matches to each other.Full size imageFigure 8The red dots (o) in each sub-figures represent the experimental data of Jin et al.1. The blue dots ((*)) are obtained from the stochastic model (3) considering: (a) The seeding 1 estimated model parameters with (D= 0.002), (Q= 0.06) and (lambda = 0.4). (b) The seeding 2 estimated model parameters with (D= 0.01), (Q= 0.15) and (lambda = 0.6). (c) The seeding 3 estimated model parameters with (D= 0.002), (Q= 0.2) and (lambda = 0.4).Full size imageFigures 7 and 8 suggest that the stochastic model is valid. So the model can be further analyzed to meet the first objective. Differentiating (P_{st} (x)), we obtain (frac{dP_{st} (x)}{dx}=frac{N^{prime }}{[B(x)]^2} exp left( int frac{A(x)}{B(x)}dx right) left( A(x)-frac{dB(x)}{dx} right)) and (frac{d^{2}P_{st} (x)}{dx^{2}}= frac{N^{prime }}{[B(x)]^{2}}exp left( int frac{A(x)}{B(x)}dx right) left( frac{dA(x)}{dx}-frac{d^{2}B(x)}{dx^{2}} right) +frac{N^{prime }}{[B(x)]^{2}} left( A(x)-frac{dB(x)}{dx} right) exp left( int frac{A(x)}{B(x)}dx right) frac{A(x)}{B(x)}-frac{2}{[B(x)]^3}N^{prime } exp left( int frac{A(x)}{B(x)}dx right) left( A(x)-frac{dB(x)}{dx} right) frac{dB(x)}{dx}). At the extrema of the SSPDF, we must have (frac{dP_{st} (x)}{dx}=0) i.e. (left( A(x)-frac{dB(x)}{dx} right) =0).

Theorem 3

(x^{*}approx K-K left( frac{nK^{delta +1}+D(delta +1) K^{2delta +1}-lambda sqrt{DQ}(delta +1)K^{delta }}{beta r K^{alpha +1}+n(delta +1) K^{(delta +1)}+D(delta +1) (2delta +1)K^{(2delta +1)}-lambda sqrt{DQ}delta (delta +1)K^{delta }} right)) is the conditional MSSCD due to the correlated additive and multiplicative noises under the condition (r_{p}(alpha +1)x^{*}{}^{alpha }-frac{r_{p}}{K^{beta }}(alpha +beta +1)x^{*}{}^{(alpha +beta )} -n(delta +1)x^{*}{}^{delta }-D(delta +1)(2delta +1)x^{*}{}^{(2delta )}+lambda sqrt{Dalpha }delta (delta +1)x^{*}{}^{(delta -1)} < 0) (proof is in the supplementary information). Figure 9 visualizes the effect of noise strength and correlation strength on the conditional MSSCD. The conditional MSSCD increases with the additive noise strength (Q) and decreases with the multiplicative noise strength (D) when the other model parameters are fixed (Fig. 9a). There is a high chance of overproliferation for a low D and a high Q (Fig. 9a). Again, there is a high chance of extinction for the low Q and high D. The conditional MSSCD depends more on D than (lambda) (Fig. 9b), and more on (lambda) than Q (Fig. 9c). The conditional MSSCD increases with (lambda) and Q; there is a high chance of overproliferation for high (lambda) and Q. The extinction risk of cells from the culture increases with low (lambda) and Q.Figure 9The change in the conditional MSSCD value for different noise strengths and correlation strength using the parameters estimated for seeding 1: (a) the conditional MSSCD values in (D-Q) noise strength plane with highest correlation ((lambda =1)); (b) the conditional MSSCD values in (D-lambda) noise plane with (Q=0.01); (c) the conditional MSSCD values in (Q-lambda) noise plane with (D=0.01).Full size imageDue to the difficulty and complicated expression of the analytical expression of the SSPDF, we use numerical simulation to study the steady-state behavior in the long run under correlated noises. We draw a histogram of the cell densities based on 500 normal sample paths at the final time points. We use seeding 1 fitting estimates as the initial parameter values for this simulation. The cell population is stable and steady at either 0 cell density or at the conditional MSSCD. The distribution is symmetric around the conditional MSSCD for (lambda =1) (Fig. 10a). There is a loss in the symmetry for the decreasing (lambda). For (lambda =0.5), there is a mode at the zero states with another mode at conditional MSSCD (Fig. 10b). The histogram shows a bi-modality for low values of (lambda). The mode at the zero state is highest for (lambda =0) (Fig. 10c). Therefore, the extinction chance increases for zero noise correlation between the additive and the multiplicative noises.Figure 10Distribution of cell density for (r_{p}=0.13), (K=1.43), (n=0.0095), (alpha =1.15), (beta =0.99), (delta =0.2), (D=0.01), (Q=0.01), and variable correlation between additive and multiplicative noises: (a) (lambda =1), (b) (lambda =0.5) and (c) (lambda =0).Full size imageThe sustainability of the cell population depends on the strength of the two noises, like the correlation strength between them. For the zero strength multiplicative noise, the population has the mode at around the conditional MSSCD value (Fig. 11). Therefore, the population sustains in this case and tries to stabilize at the conditional MSSCD value. For (D=0.02), there is a bimodality, where the highest mode is at the zero cell density. For (D=0.05), we observe only one mode at (x=0). Therefore, with the increasing values of the multiplicative noise strengths (D), the chance of extinction increases for (lambda =0.5), (Q=0.01), and other constant model parameters for the seeding condition 1. Similar things happen for increasing Q values considering (D=0.01), (lambda =0.5), and other constant model parameters (Fig. 12).Figure 11Distribution of cell density for (r_{p}=0.13), (K=1.43), (n=0.0095), (alpha =1.15), (beta =0.99), (delta =0.2), (lambda =0.5), (Q=0.01), and variable strength of multiplicative noise: (a) (D=0.05), (b) (D=0.02) and (c) (D=0).Full size imageFigure 12Distribution of cell density for (r_{p}=0.13), (K=1.43), (n=0.0095), (alpha =1.15), (beta =0.99), (delta =0.2), (lambda =0.5), (D=0.01), and variable correlation between multiplicative noise: (a) (Q=0.05), (b) (Q=0.02) and (c) (Q=0).Full size image Remark 5 We have previously discussed the scenario for (alpha =delta) for deterministic case in Remark 4. It is important to understand the scenario under stochastic case too. For (alpha =delta) the proposed stochastic model 3 becomes (frac{dx(t)}{dt}=r_{p}x(t)^{(alpha +1)}left( 1-big (frac{x(t)}{K}big )^{beta }right) - nx(t)^{(alpha +1)}-x(t)^{(alpha +1)} epsilon (t)+ Gamma (t)). For this stochastic model (g_{1}(x)=-x^{alpha +1}) and (g_{2}(x)=1). We get, (A(x)=r_{p}x^{alpha +1}left( 1-Big (frac{x}{K}Big )^{beta } right) -nx^{(alpha +1)}+D(alpha +1)x^{(2alpha +1)}-lambda sqrt{DQ}(alpha +1)x^{alpha }) and (B(x)=Dx^{(2alpha +2)}-2lambda sqrt{DQ}x^{(alpha +1)}+Q). The extrema of the SPDF (big (x(t)=x^{*}big )) must satisfy the growth equation (r_{p}{x^{*}}^{alpha +1}-frac{r_{p}}{K^{beta }}(x^{*})^{alpha +beta +1}-n(x^{*})^{alpha +1}-D(alpha +1)(x^{*})^{2alpha +1}+lambda sqrt{D~Q}(alpha +1)(x^{*})^{alpha }=0). Therefore, for (alpha =delta) the conditional MSSCD is (x^{*}= K-Kfrac{nK^{(alpha +1)}+D(alpha +1)K^{(2alpha +1)}-lambda sqrt{DQ}(alpha +1)K^{alpha }}{beta r_{p}K^{(alpha +1)}+nK^{(alpha +1)}(alpha +1)+D(alpha +1)(2alpha +1)K^{(2alpha +1)}-alpha lambda sqrt{DQ}(alpha +1)K^{alpha }}) under the condition ((r_{p}-n)(alpha +1)(x^{*})^{alpha }-frac{r_{p}}{K^{beta }}(alpha +beta +1)(x^{*})^{(alpha + beta )}-(alpha +1)(2alpha +1)D(x^{*})^{2alpha }+lambda sqrt{DQ}(alpha +1)alpha (x^{*})^{(alpha -1)} More

Speijer, D., Lukeš, J. & Eliáš, M. Sex is a ubiquitous, ancient, and inherent attribute of eukaryotic life. Proc. Natl Acad. Sci. 112, 8827–8834 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bachtrog, D. et al. Sex determination: why so many ways of doing it? PLoS Biol. 12, e1001899 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Ah-King, M. & Nylin, S. Sex in an evolutionary perspective: just another reaction norm. Evolut. Biol. 37, 234–246 (2010).Article 

Google Scholar 
Leonard, J. L. The evolution of sexual systems in animals. In: Leonard, J.L. (ed.). Transitions between sexual systems: understanding the mechanisms of, and pathways between, dioecy, hermaphroditism and other sexual systems, 1–58 Springer (2019).Weeks, S. C., Benvenuto, C. & Reed, S. K. When males and hermaphrodites coexist: a review of androdioecy in animals. Integr. Comp. Biol. 46, 449–464 (2006).PubMed 
Article 

Google Scholar 
Goldberg, E. E. et al. Macroevolutionary synthesis of flowering plant sexual systems. Evolution 71, 898–912 (2017).PubMed 
Article 

Google Scholar 
Waples, R. S., Mariani, S. & Benvenuto, C. Consequences of sex change for effective population size. Proc. R. Soc. B: Biol. Sci. 285, 20181702 (2018).Article 

Google Scholar 
Benvenuto, C. & Weeks, S. C. Hermaphroditism and gonochorism. The Natural History of the Crustacea: Reproductive Biology VI, 197–241 (2020).
Google Scholar 
Mariani, S., Sala-Bozano, M., Chopelet, J. & Benvenuto, C. Spatial and temporal patterns of size-at-sex-change in two exploited coastal fish. Environ. Biol. Fishes 96, 535–541 (2013).Article 

Google Scholar 
Käfer, J., Marais, G. A. & Pannell, J. R. On the rarity of dioecy in flowering plants. Mol. Ecol. 26, 1225–1241 (2017).PubMed 
Article 

Google Scholar 
Atz, J. Intersexuality in Fishes. In C.N. Amstrong and A.J. Marshall (eds). Intersexuality in vertebrates including man, 145–232 Academic Press, London (1964).Jarne, P. & Auld, J. R. Animals mix it up too: the distribution of self-fertilization among hermaphroditic animals. Evolution 60, 1816–1824 (2006).PubMed 
Article 

Google Scholar 
Leonard, J. L. Williams’ paradox and the role of phenotypic plasticity in sexual systems. Integr. Comp. Biol. 53, 671–688 (2013).PubMed 
Article 

Google Scholar 
Weeks, S. C. The role of androdioecy and gynodioecy in mediating evolutionary transitions between dioecy and hermaphroditism in the animalia. Evolution 66, 3670–3686 (2012).PubMed 
Article 

Google Scholar 
Renner, S. S. The relative and absolute frequencies of angiosperm sexual systems: dioecy, monoecy, gynodioecy, and an updated online database. Am. J. Bot. 101, 1588–1596 (2014).PubMed 
Article 

Google Scholar 
Bawa, K. S. Evolution of dioecy in flowering plants. Annu. Rev. Ecol. Syst. 11, 15–39 (1980).Article 

Google Scholar 
Charlesworth, B. & Charlesworth, D. A model for the evolution of dioecy and gynodioecy. Am. Nat. 112, 975–997 (1978).Article 

Google Scholar 
Charlesworth, D. Androdioecy and the evolution of dioecy. Biol. J. Linn. Soc. 22, 333–348 (1984).Article 

Google Scholar 
Pannell, J. R. The evolution and maintenance of androdioecy. In: Annual Review of Ecology and Systematics 397–425 (2002).Bull, J. & Charnov, E. On irreversible evolution. Evolution 39, 1149–1155 (1985).CAS 
PubMed 
Article 

Google Scholar 
Barrett, S. C. The evolution of plant reproductive systems: how often are transitions irreversible? Proc. R. Soc. B: Biol. Sci. 280, 20130913 (2013).Article 

Google Scholar 
Oyarzún, P. A., Nuñez, J. J., Toro, J. E. & Gardner, J. P. Trioecy in the marine mussel Semimytilus algosus (Mollusca, Bivalvia): stable sex ratios across 22 degrees of a latitudinal gradient. Front. Mar. Sci. 7, 348 (2020).Article 

Google Scholar 
Dani, K. & Kodandaramaiah, U. Plant and animal reproductive strategies: lessons from offspring size and number tradeoffs. Front. Ecol. Evol. 5, 38 (2017).Article 

Google Scholar 
Avise, J. & Mank, J. Evolutionary perspectives on hermaphroditism in fishes. Sex. Dev. 3, 152–163 (2009).CAS 
PubMed 
Article 

Google Scholar 
Dornburg, A. & Near, T. J. The Emerging phylogenetic perspective on the evolution of Actinopterygian fishes. Annu. Rev. Ecol. Evol. Syst. 52, 427–452 (2021).Article 

Google Scholar 
Costa, W. J., Lima, S. M. & Bartolette, R. Androdioecy in Kryptolebias killifish and the evolution of self-fertilizing hermaphroditism. Biol. J. Linn. Soc. 99, 344–349 (2010).Article 

Google Scholar 
Costa, W. Colouration, taxonomy and geographical distribution of mangrove killifishes, the Kryptolebias marmoratus species group, in southern Atlantic coastal plains of Brazil (Cyprinodontiformes: Rivulidae). Ichthyol. Explor. Freshw. 27, 183–192 (2016).
Google Scholar 
Powell, M. L., Kavanaugh, S. I. & Sower, S. A. Seasonal concentrations of reproductive steroids in the gonads of the Atlantic hagfish, Myxine glutinosa. J. Exp. Zool. Part A Comp. Exp. Biol. 301, 352–360 (2004).Article 
CAS 

Google Scholar 
Pennell, M. W., Mank, J. E. & Peichel, C. L. Transitions in sex determination and sex chromosomes across vertebrate species. Mol. Ecol. 27, 3950–3963 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Ghiselin, M. T. The evolution of hermaphroditism among animals. Q. Rev. Biol. 44, 189–208 (1969).CAS 
PubMed 
Article 

Google Scholar 
Eppley, S. M. & Jesson, L. K. Moving to mate: the evolution of separate and combined sexes in multicellular organisms. J. Evol. Biol. 21, 727–736 (2008).CAS 
PubMed 
Article 

Google Scholar 
Warner, R. R. The adaptive significance of sequential hermaphroditism in animals. Am. Nat. 109, 61–82 (1975).Article 

Google Scholar 
Warner, R. R., Robertson, D. R. & Leigh, E. G. Sex change and sexual selection. Science 190, 633–638 (1975).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Charnov, E. L. The Theory of Sex Allocation. Princeton University Press, USA (1982).Policansky, D. Sex change in plants and animals. Annu. Rev. Ecol. Syst. 13, 471–495 (1982).Article 

Google Scholar 
Benvenuto, C., Coscia, I., Chopelet, J., Sala-Bozano, M. & Mariani, S. Ecological and evolutionary consequences of alternative sex-change pathways in fish. Sci. Rep. 7, 9084 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Charnov, E. L. Natural selection and sex change in pandalid shrimp: test of a life-history theory. Am. Nat. 113, 715–734 (1979).MathSciNet 
Article 

Google Scholar 
Broquet, T. et al. The size advantage model of sex allocation in the protandrous sex-changer Crepidula fornicata: role of the mating system, sperm storage, and male mobility. Am. Nat. 186, 404–420 (2015).PubMed 
Article 

Google Scholar 
Erisman, B. E., Craig, M. T. & Hastings, P. A. A phylogenetic test of the size-advantage model: evolutionary changes in mating behavior influence the loss of sex change in a fish lineage. Am. Nat. 174, E83–E99 (2009).PubMed 
Article 

Google Scholar 
Buxton, C. D. & Garratt, P. A. Alternative reproductive styles in seabreams (Pisces: Sparidae). Environ. Biol. Fishes 28, 113–124 (1990).Article 

Google Scholar 
Shapiro, D. Y. Social behavior, group structure, and the control of sex reversal in hermaphroditic fish. Adv. Study Behav. 10, 43–102 (1979).Article 

Google Scholar 
Stearns, S. C. Life history evolution: successes, limitations, and prospects. Naturwissenschaften 87, 476–486 (2000).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Waples, R. S., Luikart, G., Faulkner, J. R. & Tallmon, D. A. Simple life-history traits explain key effective population size ratios across diverse taxa. Proc. R. Soc. Lond. B: Biol. Sci. 280, 20131339 (2013).
Google Scholar 
Martinez, A. S., Willoughby, J. R. & Christie, M. R. Genetic diversity in fishes is influenced by habitat type and life-history variation. Ecol. Evolution 8, 12022–12031 (2018).Article 

Google Scholar 
Harvey, P. H. & Pagel, M. D. The comparative method in evolutionary biology. (Oxford University Press, USA, 1991).Barneche, D. R., Robertson, D. R., White, C. R. & Marshall, D. J. Fish reproductive-energy output increases disproportionately with body size. Science 360, 642–645 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Brandl, S. J. & Bellwood, D. R. Pair-formation in coral reef fishes: an ecological perspective. Oceanogr. Mar. Biol.: Annu. Rev. 52, 1–80 (2014).
Google Scholar 
Fitzpatrick, J. L. Sperm competition and fertilization mode in fishes. Philos. Trans. R. Soc. B: Biol. Sci. 375, 20200074 (2020).Article 

Google Scholar 
Parker, G. A. Conceptual developments in sperm competition: a very brief synopsis. Philos. Trans. R. Soc. B: Biol. Sci. 375, 20200061 (2020).Article 

Google Scholar 
Warner, R. R. Sex change in fishes: hypotheses, evidence, and objections. Environ. Biol. Fishes 22, 81–90 (1988).Article 

Google Scholar 
Molloy, P. P., Goodwin, N. B., Côté, I. M., Reynolds, J. D. & Gage, M. J. Sperm competition and sex change: a comparative analysis across fishes. Evolution 61, 640–652 (2007).PubMed 
Article 

Google Scholar 
Erisman, B. E., Petersen, C. W., Hastings, P. A. & Warner, R. R. Phylogenetic perspectives on the evolution of functional hermaphroditism in teleost fishes. Integr. Comp. Biol. 53, 736–754 (2013).PubMed 
Article 

Google Scholar 
Sadovy, Y., Colin, P. & Domeier, M. Aggregation and spawning in the tiger grouper, Mycteroperca tigris (Pisces: Serranidae). Copeia 1994, 511–516 (1994).Article 

Google Scholar 
Muñoz, R. C. & Warner, R. R. A new version of the size-advantage hypothesis for sex change: incorporating sperm competition and size-fecundity skew. Am. Nat. 161, 749–761 (2003).PubMed 
Article 

Google Scholar 
Horne, C. R., Hirst, A. G. & Atkinson, D. Selection for increased male size predicts variation in sexual size dimorphism among fish species. Proc. R. Soc. B: Biol. Sci. 287, 20192640 (2020).Article 

Google Scholar 
Parker, G. The evolution of expenditure on testes. J. Zool. 298, 3–19 (2016).Article 

Google Scholar 
Stockley, P., Gage, M., Parker, G. & Møller, A. Sperm competition in fishes: the evolution of testis size and ejaculate characteristics. Am. Nat. 149, 933–954 (1997).CAS 
PubMed 
Article 

Google Scholar 
Pla, S., Benvenuto, C., Capellini, I. & Piferrer, F. A phylogenetic comparative analysis on the evolution of sequential hermaphroditism in seabreams (Teleostei: Sparidae). Sci. Rep. 10, 3606 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Vrijenhoek, R. C. Unisexual fish: model systems for studying ecology and evolution. Annu. Rev. Ecol. Syst. 25, 71–96 (1994).Article 

Google Scholar 
Sadovy de Mitcheson, Y. & Liu, M. Functional hermaphroditism in teleosts. Fish. Fish. 9, 1–43 (2008).Article 

Google Scholar 
Rabosky, D. L. et al. An inverse latitudinal gradient in speciation rate for marine fishes. Nature 559, 392 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Froese, R., Pauly, D. & Editors. FishBase. World Wide Web electronic publication. www.fishbase.org (2018).Moore, W. S. Evolutionary ecology of unisexual fishes. In: Evolutionary genetics of fishes, 329–398 (Springer, 1984).Collin, R. & Miglietta, M. P. Reversing opinions on Dollo’s Law. Trends Ecol. Evol. 23, 602–609 (2008).PubMed 
Article 

Google Scholar 
Domes, K., Norton, R. A., Maraun, M. & Scheu, S. Re-evolution of sexuality breaks Dollo’s law. Proc. Natl Acad. Sci. 104, 7139–7144 (2007).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Dollo, L. Les lois de l’évolution. Bull. Soc. Belge Géol. Paléont. Hydrol. 7, 164–166 (1893).
Google Scholar 
King, B. & Lee, M. S. Ancestral state reconstruction, rate heterogeneity, and the evolution of reptile viviparity. Syst. Biol. 64, 532–544 (2015).CAS 
PubMed 
Article 

Google Scholar 
Uller, T. & Helanterä, H. From the origin of sex-determining factors to the evolution of sex-determining systems. Q. Rev. Biol. 86, 163–180 (2011).PubMed 
Article 

Google Scholar 
Devlin, R. H. & Nagahama, Y. Sex determination and sex differentiation in fish: an overview of genetic, physiological, and environmental influences. Aquaculture 208, 191–364 (2002).CAS 
Article 

Google Scholar 
Volff, J.-N., Nanda, I., Schmid, M. & Schartl, M. Governing sex determination in fish: regulatory putsches and ephemeral dictators. Sex. Dev. 1, 85–99 (2007).PubMed 
Article 

Google Scholar 
Nagahama, Y., Chakraborty, T., Paul-Prasanth, B., Ohta, K. & Nakamura, M. Sex determination, gonadal sex differentiation and plasticity in vertebrate species. Physiol. Rev. 101, 1237–1308 (2020).PubMed 
Article 

Google Scholar 
Penman, D. J. & Piferrer, F. Fish gonadogenesis. Part I: genetic and environmental mechanisms of sex determination. Rev. Fish. Sci. 16(S1), 16–34 (2008).CAS 
Article 

Google Scholar 
Mank, J. E., Promislow, D. E. L. & Avise, J. C. Evolution of alternative sex-determining mechanisms in teleost fishes. Biol. J. Linn. Soc. 87, 83–93 (2006).Article 

Google Scholar 
Galetti, P. M., Aguilar, C. T. & Molina, W. F. An overview of marine fish cytogenetics. Hydrobiologia 420, 55–62 (2000).Article 

Google Scholar 
Yoshida, K. et al. Sex chromosome turnover contributes to genomic divergence between incipient stickleback species. PLoS Genet. 10, e1004223 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Ross, J. A., Urton, J. R., Boland, J., Shapiro, M. D. & Peichel, C. L. Turnover of sex chromosomes in the stickleback fishes (Gasterosteidae). PLoS Genet. 5, e1000391 (2009).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Vicoso, B. Molecular and evolutionary dynamics of animal sex-chromosome turnover. Nature Ecology & Evolution 1–10 (2019).Gamble, T. et al. Restriction site-associated DNA sequencing (RAD-seq) reveals an extraordinary number of transitions among gecko sex-determining systems. Mol. Biol. Evol. 32, 1296–1309 (2015).CAS 
PubMed 
Article 

Google Scholar 
Pokorná, M. & Kratochvíl, L. Phylogeny of sex-determining mechanisms in squamate reptiles: are sex chromosomes an evolutionary trap? Zool. J. Linn. Soc. 156, 168–183 (2009).Article 

Google Scholar 
Furman, B. L. et al. Sex chromosome evolution: sso many exceptions to the rules. Genome Biol. Evol. 12, 750–763 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Carvalho, N. D. M. et al. Cytogenetics of Synbranchiformes: a comparative analysis of two Synbranchus Bloch, 1795 species from the Amazon. Genetica 140, 149–158 (2012).CAS 
PubMed 
Article 

Google Scholar 
Piferrer, F. Epigenetic mechanisms in sex determination and in the evolutionary transitions between sexual systems. Philos. Trans. R. Soc. B: Biol. Sci. 376, 20200110 (2021).Article 
CAS 

Google Scholar 
Grant, S. et al. Genetics of sex determination in flowering plants. Dev. Genet. 15, 214–230 (1994).Article 

Google Scholar 
Harrington Jr, R. W. How ecological and genetic factors interact to determine when self-fertilizing hermaphrodites of Rivulus marmoratus change into functional secondary males, with a reappraisal of the modes of intersexuality among fishes. Copeia 389–432 (1971).Adolfi, M. C., Nakajima, R. T., Nóbrega, R. H. & Schartl, M. Intersex, Hermaphroditism, and gonadal plasticity in vertebrates: Evolution of the Müllerian duct and Amh/Amhr2 signalling. Annual Review of Animal Biosciences (2018).Adkins-Regan, E. Early organizational effects of hormones: an evolutionary perspective. In Adler, N.T. (ed.) Neuroendocrinology of reproduction: physiology and behavior, 159–228 (Springer, USA, 1981).Navara, K. J. The truth about Nemo’s dad: sex-changing behaviors in fishes. In Choosing Sexes 183–212 (Springer, Cham, 2018).Orban, L., Sreenivasan, R. & Olsson, P. E. Long and winding roads: testis differentiation in zebrafish. Mol. Cell. Endocrinol. 312, 35–41 (2009).CAS 
PubMed 
Article 

Google Scholar 
Zohar, Y., Abraham, M. & Gordin, H. The gonadal cycle of the captivity-reared hermaphroditic teleost Sparus aurata (L.) during the first two years of life. Annales de. Biologie Anim. Biochim. Biophys. 18, 877–882 (1978).Article 

Google Scholar 
Chang, C.-F. & Yueh, W.-S. Annual cycle of gonadal histology and steroid profiles in the juvenile males and adult females of the protandrous black porgy, Acanthopagrus schlegelii. Aquaculture 91, 179–196 (1990).CAS 
Article 

Google Scholar 
Miura, S., Nakamura, S., Kobayashi, Y., Piferrer, F. & Nakamura, M. Differentiation of ambisexual gonads and immunohistochemical localization of P450 cholesterol side-chain cleavage enzyme during gonadal sex differentiation in the protandrous anemonefish, Amphiprion clarkii. Comp. Biochem. Physiol. Part B: Biochem. Mol. Biol. 149, 29–37 (2008).Article 
CAS 

Google Scholar 
Yamaguchi, S. & Iwasa, Y. Advantage for the sex changer who retains the gonad of the nonfunctional sex. Behav. Ecol. Sociobiol. 71, 39 (2017).Article 

Google Scholar 
Munday, P. L., Kuwamura, T. & Kroon, F. J. Bi-directional sex change in marine fishes. In: Cole, K.S. (ed.) Reproduction and sexuality in marine fishes: Patterns and processes. 241–271 (University of California Press, Berkeley, USA, 2010).Uller, T., Feiner, N., Radersma, R., Jackson, I. S. & Rago, A. Developmental plasticity and evolutionary explanations. Evol. Dev. 22, 47–55 (2020).PubMed 
Article 

Google Scholar 
Pla, S., Maynou, F. & Piferrer, F. Hermaphroditism in fish: incidence, distribution and associations with abiotic environmental factors. Rev. Fish. Biol. Fish. 31, 935–955 (2021).Article 

Google Scholar 
Boettiger, C., Lang, D. T. & Wainwright, P. C. rfishbase: exploring, manipulating and visualizing FishBase data from R. J. Fish. Biol. 81, 2030–2039 (2012).CAS 
PubMed 
Article 

Google Scholar 
Pagel, M., Meade, A. & Barker, D. Bayesian estimation of ancestral character states on phylogenies. Syst. Biol. 53, 673–684 (2004).PubMed 
Article 

Google Scholar 
Pagel, M. & Meade, A. Bayesian analysis of correlated evolution of discrete characters by reversible-jump Markov chain Monte Carlo. Am. Nat. 167, 808–825 (2006).PubMed 
Article 

Google Scholar 
Pagel, M. Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete. Proc. R. Soc. B: Biol. Sci. 255, 37–45 (1994).ADS 
Article 

Google Scholar 
Currie, T. E. & Meade, A. In Modern phylogenetic comparative methods and their application in evolutionary biology, 263–286 (Springer, 2014).Furness, A. I. & Capellini, I. The evolution of parental care diversity in amphibians. Nat. Commun. 10, 1–12 (2019).CAS 
Article 

Google Scholar 
Rambaut, A., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67, 901 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
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 
Freckleton, R. P., Harvey, P. H. & Pagel, M. Phylogenetic analysis and comparative data: a test and review of evidence. Am. Nat. 160, 712–726 (2002).CAS 
PubMed 
Article 

Google Scholar 
Pagel, M. Inferring evolutionary processes from phylogenies. Zool. Scr. 26, 331–348 (1997).Article 

Google Scholar 
Pagel, M. Inferring the historical patterns of biological evolution. Nature 401, 877–884 (1999).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Orme, D. The caper package: comparative analysis of phylogenetics and evolution in R. https://cran.r-project.org/web/packages/caper/vignettes/caper.pdf (2018).Schiettekatte, N., Brandl, S. & Casey, J. Fishualize: Color palettes based on fish species. R package v0.2.2 (2021). More

Flemming, H. C. & Wuertz, S. Bacteria and archaea on Earth and their abundance in biofilms. Nat. Rev. Microbiol. 17, 247–260 (2019).CAS 
PubMed 
Article 

Google Scholar 
Bulgarelli, D., Schlaeppi, K., Spaepen, S., Ver Loren van Themaat, E. & Schulze-Lefert, P. Structure and functions of the bacterial microbiota of plants. Annu Rev. Plant Biol. 64, 807–838 (2013).CAS 
PubMed 
Article 

Google Scholar 
Turnbaugh, P. J. et al. The human microbiome project. Nature 449, 804–810 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Venturelli, O. S. et al. Deciphering microbial interactions in synthetic human gut microbiome communities. Mol. Syst. Biol. 14, e8157 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Foster, K. R. & Bell, T. Competition, not cooperation, dominates interactions among culturable microbial species. Curr. Biol. 22, 1845–1850 (2012).CAS 
PubMed 
Article 

Google Scholar 
Helfrich, E. J. N. et al. Bipartite interactions, antibiotic production and biosynthetic potential of the Arabidopsis leaf microbiome. Nat. Microbiol. 3, 909–919 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Coyte, K. Z. & Rakoff-Nahoum, S. Understanding competition and cooperation within the mammalian gut microbiome. Curr. Biol. 29, 538–544 (2019).Article 
CAS 

Google Scholar 
Turner, T. R. et al. Comparative metatranscriptomics reveals kingdom level changes in the rhizosphere microbiome of plants. ISME J. 7, 2248–2258 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Trivedi, P., Leach, J. E., Tringe, S. G., Sa, T. & Singh, B. K. Plant–microbiome interactions: from community assembly to plant health. Nat. Rev. Microbiol. 18, 607–621 (2020).CAS 
PubMed 
Article 

Google Scholar 
Müller, D. B., Vogel, C., Bai, Y. & Vorholt, J. A. The plant microbiota: systems-level insights and perspectives. Annu. Rev. Genet. 50, 211–234 (2016).PubMed 
Article 
CAS 

Google Scholar 
Lugtenberg, B. & Kamilova, F. Plant-growth-promoting Rhizobacteria. Annu. Rev. Microbiol. 63, 541–556 (2009).CAS 
PubMed 
Article 

Google Scholar 
Berendsen, R. L., Pieterse, C. M. J. & Bakker, P. A. H. M. The rhizosphere microbiome and plant health. Trends Plant Sci. 17, 478–486 (2012).CAS 
PubMed 
Article 

Google Scholar 
Innerebner, G., Knief, C. & Vorholt, J. A. Protection of Arabidopsis thaliana against leaf-pathogenic Pseudomonas syringae by Sphingomonas strains in a controlled model system. Appl. Environ. Microbiol. 77, 3202–3210 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Shekhawat, K. et al. Root endophyte induced plant thermotolerance by constitutive chromatin modification at heat stress memory gene loci. EMBO Rep. 22, e51049 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Vorholt, J. A. Microbial life in the phyllosphere. Nat. Rev. Microbiol. 10, 828–840 (2012).CAS 
PubMed 
Article 

Google Scholar 
Bodenhausen, N., Bortfeld-Miller, M., Ackermann, M. & Vorholt, J. A. A synthetic community approach reveals plant genotypes affecting the phyllosphere microbiota. PLoS Genet. 10, e1004283 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Reisberg, E. E., Hildebrandt, U., Riederer, M. & Hentschel, U. Distinct phyllosphere bacterial communities on Arabidopsis wax mutant leaves. PLoS ONE 8, e78613 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kniskern, J. M., Traw, M. B. & Bergelson, J. Salicylic acid and jasmonic acid signaling defense pathways reduce natural bacterial diversity on Arabidopsis thaliana. Mol. Plant Microbe Interact. 20, 1512–1522 (2007).CAS 
PubMed 
Article 

Google Scholar 
Pfeilmeier, S. et al. The plant NADPH oxidase RBOHD is required for microbiota homeostasis in leaves. Nat. Microbiol. 6, 852–864 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Chen, T. et al. A plant genetic network for preventing dysbiosis in the phyllosphere. Nature 580, 653–657 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Hassani, M. A., Duran, P. & Hacquard, S. Microbial interactions within the plant holobiont. Microbiome 6, 58 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Lidicker, W. Z. Clarification of interactions in ecological systems. Bioscience 29, 475–477 (1979).Article 

Google Scholar 
Schlechter, R. O., Miebach, M. & Remus-Emsermann, M. N. P. Driving factors of epiphytic bacterial communities: a review. J. Adv. Res. 19, 57–65 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Faust, K. & Raes, J. Microbial interactions: from networks to models. Nat. Rev. Microbiol. 10, 538–550 (2012).CAS 
PubMed 
Article 

Google Scholar 
Grosskopf, T. & Soyer, O. S. Synthetic microbial communities. Curr. Opin. Microbiol. 18, 72–77 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Blair, P. M. et al. Exploration of the biosynthetic potential of the Populus microbiome. mSystems 3, e00045-00018 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Suda, W., Nagasaki, A. & Shishido, M. Powdery mildew-infection changes bacterial community composition in the phyllosphere. Microbes Environ. 24, 217–223 (2009).PubMed 
Article 

Google Scholar 
Manching, H. C., Balint-Kurti, P. J. & Stapleton, A. E. Southern leaf blight disease severity is correlated with decreased maize leaf epiphytic bacterial species richness and the phyllosphere bacterial diversity decline is enhanced by nitrogen fertilization. Front. Plant Sci. 5, 403 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Agler, M. T. et al. Microbial hub taxa link host and abiotic factors to plant microbiome variation. PLoS Biol. 14, 100235 (2016).Article 
CAS 

Google Scholar 
Layeghifard, M., Hwang, D. M. & Guttman, D. S. Disentangling interactions in the microbiome: a network perspective. Trends Microbiol. 25, 217–228 (2017).CAS 
PubMed 
Article 

Google Scholar 
Faust, K. et al. Microbial co-occurrence relationships in the human microbiome. PLoS Comput. Biol. 8, e1002606 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Carr, A., Diener, C., Baliga, N. S. & Gibbons, S. M. Use and abuse of correlation analyses in microbial ecology. ISME J. 13, 2647–2655 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Vorholt, J. A., Vogel, C., Carlström, C. I. & Müller, D. B. Establishing causality: opportunities of synthetic communities for plant microbiome research. Cell Host Microbe 22, 142–155 (2017).CAS 
PubMed 
Article 

Google Scholar 
Bai, Y. et al. Functional overlap of the Arabidopsis leaf and root microbiota. Nature 528, 364–369 (2015).CAS 
PubMed 
Article 

Google Scholar 
Knief, C., Frances, L. & Vorholt, J. A. Competitiveness of diverse Methylobacterium strains in the phyllosphere of Arabidopsis thaliana and identification of representative models, including M. extorquens PA1. Microb. Ecol. 60, 440–452 (2010).PubMed 
Article 

Google Scholar 
Fan, J., Crooks, C. & Lamb, C. High-throughput quantitative luminescence assay of the growth in planta of Pseudomonas syringae chromosomally tagged with Photorhabdus luminescens luxCDABE. Plant J. 53, 393–399 (2008).CAS 
PubMed 
Article 

Google Scholar 
Carlström, C. I. et al. Synthetic microbiota reveal priority effects and keystone strains in the Arabidopsis phyllosphere. Nat. Ecol. Evol. 3, 1445–1454 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Vogel, C. M., Potthoff, D. M., Schäfer, M., Barandun, N. & Vorholt, J. A. Protective role of the Arabidopsis leaf microbiota against a bacterial pathogen. Nat. Microbiol. 6, 1537–1548 (2021).CAS 
PubMed 
Article 

Google Scholar 
Chen, I.-M. A. et al. The IMG/M data management and analysis system v.6.0: new tools and advanced capabilities. Nucleic Acids Res. 49, 751–763 (2020).Article 
CAS 

Google Scholar 
Ortiz, A., Vega, N. M., Ratzke, C. & Gore, J. Interspecies bacterial competition regulates community assembly in the C. elegans intestine. ISME J. 15, 2131–2145 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Goberna, M. & Verdú, M. Predicting microbial traits with phylogenies. ISME J. 10, 959–967 (2016).PubMed 
Article 

Google Scholar 
Webb, C. O., Ackerly, D. D., McPeek, M. A. & Donoghue, M. J. Phylogenies and community ecology. Annu. Rev. Ecol. Syst. 33, 475–505 (2002).Article 

Google Scholar 
Cahill, J. F., Kembel, S. W., Lamb, E. G. & Keddy, P. A. Does phylogenetic relatedness influence the strength of competition among vascular plants? Perspect. Plant Ecol. 10, 41–50 (2008).Article 

Google Scholar 
Maherali, H. & Klironomos, J. N. Influence of phylogeny on fungal community assembly and ecosystem functioning. Science 316, 1746–1748 (2007).CAS 
PubMed 
Article 

Google Scholar 
Duncan, R. P. & Williams, P. A. Ecology – Darwin’s naturalization hypothesis challenged. Nature 417, 608–609 (2002).CAS 
PubMed 
Article 

Google Scholar 
Slingsby, J. A. & Verboom, G. A. Phylogenetic relatedness limits co-occurrence at fine spatial scales: evidence from the schoenoid sedges (Cyperaceae: Schoeneae) of the Cape Floristic Region, South Africa. Am. Nat. 168, 14–27 (2006).PubMed 
Article 

Google Scholar 
Mayfield, M. M. & Levine, J. M. Opposing effects of competitive exclusion on the phylogenetic structure of communities. Ecol. Lett. 13, 1085–1093 (2010).PubMed 
Article 

Google Scholar 
Teixeira, P. J. P. L., Colaianni, N. R., Fitzpatrick, C. R. & Dangl, J. L. Beyond pathogens: microbiota interactions with the plant immune system. Curr. Opin. Microbiol. 49, 7–17 (2019).CAS 
PubMed 
Article 

Google Scholar 
Maier, B. A. et al. A general non-self response as part of plant immunity. Nat. Plants 7, 696–705 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Friedman, J., Higgins, L. M. & Gore, J. Community structure follows simple assembly rules in microbial microcosms. Nat. Ecol. Evol. 1, 0109 (2017).Article 

Google Scholar 
Kehe, J. et al. Positive interactions are common among culturable bacteria. Sci. Adv. 7, eabi7159 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Lindow, S. E. & Brandl, M. T. Microbiology of the phyllosphere. Appl. Environ. Microbiol. 69, 1875–1883 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Remus-Emsermann, M. N. P. et al. Spatial distribution analyses of natural phyllosphere-colonizing bacteria on Arabidopsis thaliana revealed by fluorescence in situ hybridization. Environ. Microbiol. 16, 2329–2340 (2014).CAS 
PubMed 
Article 

Google Scholar 
Billick, I. & Case, T. J. Higher-order interactions in ecological communities – what are they and how can they be detected. Ecology 75, 1529–1543 (1994).Article 

Google Scholar 
Grilli, J., Barabas, G., Michalska-Smith, M. J. & Allesina, S. Higher-order interactions stabilize dynamics in competitive network models. Nature 548, 210–213 (2017).CAS 
PubMed 
Article 

Google Scholar 
Levine, J. M., Bascompte, J., Adler, P. B. & Allesina, S. Beyond pairwise mechanisms of species coexistence in complex communities. Nature 546, 56–64 (2017).CAS 
PubMed 
Article 

Google Scholar 
Sundarraman, D. et al. Higher-order interactions dampen pairwise competition in the zebrafish gut microbiome. mBio 11, e01667-20 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Morris, C. in Encyclopedia for Life Sciences (National Publishing Group, 2002).Raaijmakers, J. M. & Mazzola, M. Diversity and natural functions of antibiotics produced by beneficial and plant pathogenic bacteria. Annu. Rev. Phytopathol. 50, 403–424 (2012).CAS 
PubMed 
Article 

Google Scholar 
Iversen, O. J. & Grov, A. Studies on lysostaphin – separation and characterization of 3 enzymes. Eur. J. Biochem. 38, 293–300 (1973).CAS 
PubMed 
Article 

Google Scholar 
Recsei, P. A., Gruss, A. D. & Novick, R. P. Cloning, sequence, and expression of the lysostaphin gene from Staphylococcus simulans. Proc. Natl Acad. Sci. USA 84, 1127–1131 (1987).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kessler, E., Safrin, M., Abrams, W. R., Rosenbloom, J. & Ohman, D. E. Inhibitors and specificity of Pseudomonas aeruginosa LasA. J. Biol. Chem. 272, 9884–9889 (1997).CAS 
PubMed 
Article 

Google Scholar 
Trayer, H. R. & Buckley, C. E. Molecular properties of lysostaphin, a bacteriolytic agent specific for Staphylococcus aureus. J. Biol. Chem. 245, 4842–4846 (1970).CAS 
PubMed 
Article 

Google Scholar 
Heymer, B. & Schmidt, W. C. Purification and characterization of a Streptomyces albus endo-N-acetylmuramidase lytic for group A and other beta hemolytic streptococci. Microbios 12, 51–66 (1975).CAS 
PubMed 

Google Scholar 
Vollmer, W., Joris, B., Charlier, P. & Foster, S. Bacterial peptidoglycan (murein) hydrolases. FEMS Microbiol. Rev. 32, 259–286 (2008).CAS 
PubMed 
Article 

Google Scholar 
Peyraud, R. et al. Demonstration of the ethylmalonyl-CoA pathway by using C-13 metabolomics. Proc. Natl Acad. Sci. USA 106, 4846–4851 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Schlesier, B., Breton, F. & Mock, H. P. A hydroponic culture system for growing Arabidopsis thaliana plantlets under sterile conditions. Plant Mol. Biol. Rep. 21, 449–456 (2003).CAS 
Article 

Google Scholar 
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 
Revell, L. J. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3, 217–223 (2012).Article 

Google Scholar 
Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26, 2460–2461 (2010).CAS 
PubMed 
Article 

Google Scholar 
Pruesse, E. et al. SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids Res. 35, 7188–7196 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Integrated Development Environment for R (R Studio, 2020).R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2021).Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Oksanen, J. et al. vegan: Community Ecology Package. R package v. 2.5-7 (2020).Armenteros, J. J. A. et al. SignalP 5.0 improves signal peptide predictions using deep neural networks. Nat. Biotechnol. 37, 420–423 (2019).Article 
CAS 

Google Scholar 
Gasteiger, E. et al. in The Proteomics Protocols Handbook 571–607 (ed Walker, J. M.) (Humana Press, 2005).Schindelin, J. et al. Fiji: an open-source platform for biological-image analysis. Nat. Methods 9, 676–682 (2012).CAS 
PubMed 
Article 

Google Scholar 
Bushnell, B. BBMap short read aligner, and other bioinformatic tools (SourceForge, version 38.87); https://sourceforge.net/projects/bbmapDeatherage, D. E. & Barrick, J. E. Identification of mutations in laboratory-evolved microbes from next-generation sequencing data using breseq. Methods Mol. Biol. 1151, 165–188 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Kolmogorov, M., Yuan, J., Lin, Y. & Pevzner, P. A. Assembly of long, error-prone reads using repeat graphs. Nat. Biotechnol. 37, 540–546 (2019).CAS 
PubMed 
Article 

Google Scholar 
Walker, B. J. et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS ONE 9, e112963 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Seemann, T. Prokka: rapid prokaryotic genome annotation. Bioinformatics 30, 2068–2069 (2014).CAS 
PubMed 
Article 

Google Scholar  More

Rendell, L. & Whitehead, H. Culture in whales and dolphins. Behav. Brain Sci. 24, 309–324 (2001).CAS 
Article 

Google Scholar 
Krützen, M. et al. Cultural transmission of tool use in bottlenose dolphins. Proc. Natl. Acad. Sci. U.S.A. 102, 8939–8943 (2005).ADS 
Article 

Google Scholar 
Kawai, M. Newly-acquired pre-cultural behavior of the natural troop of Japanese monkeys on Koshima Islet. Primates 6, 1–30 (1965).Article 

Google Scholar 
Slater, P. The cultural transmission of bird song. Trends Ecol. Evol. 1, 94–97 (1986).CAS 
Article 

Google Scholar 
Whitehead, H. & Rendell, L. The cultural lives of whales and dolphins. (University of Chicago Press, 2014).Whiten, A. The identification and differentiation of culture in chimpanzees and other animals: from natural history to diffusion experiments. The question of animal culture, 99–124 (2009).Allen, J. A. Community through culture: from insects to whales: How social learning and culture manifest across diverse animal communities. BioEssays 41, 1900060 (2019).Article 

Google Scholar 
Allen, J., Weinrich, M., Hoppitt, W. & Rendell, L. Network-based diffusion analysis reveals cultural transmission of lobtail feeding in humpback whales. Science 340, 485–488 (2013).ADS 
CAS 
Article 

Google Scholar 
Baker, C. Migratory movement and population structure of humpback whales (Megaptera novaeangliae) in the central and eastern North Pacific. Mar Ecol Prog Ser 31, 105–119 (1986).ADS 
Article 

Google Scholar 
Garrigue, C. et al. Movement of individual humpback whales between wintering grounds of Oceania (South Pacific), 1999 to 2004. J. Cetacean Res. Manage 3, 275–281 (2011).
Google Scholar 
Rosenbaum, H. C. et al. First circumglobal assessment of Southern Hemisphere humpback whale mitochondrial genetic variation and implications for management. Endang. Spec. Res. 32, 551–567 (2017).Article 

Google Scholar 
Noad, M. J., Cato, D. H., Bryden, M. M., Jenner, M. N. & Jenner, K. C. Cultural revolution in whale songs. Nature 408, 537. https://doi.org/10.1038/35046199 (2000).ADS 
CAS 
Article 
PubMed 

Google Scholar 
Payne, R. S. & McVay, S. Songs of humpback whales. Science 173, 585–597 (1971).ADS 
CAS 
Article 

Google Scholar 
Garland, E. C. et al. Dynamic horizontal cultural transmission of humpback whale song at the ocean basin scale. Curr. Biol. 21, 687–691. https://doi.org/10.1016/j.cub.2011.03.019 (2011).CAS 
Article 
PubMed 

Google Scholar 
Payne, R. & Guinee, L. N. Humpback whale (Megaptera novaeangliae) songs as an indicator of “stocks”. Communication and behavior of whales, 333–358 (1983).Garrigue, C. et al. Movements of humpback whales in Oceania, South Pacific. J. Cetac. Res. Manage. 4, 255–260 (2002).
Google Scholar 
Derville, S., Torres, L. G., Zerbini, A. N., Oremus, M. & Garrigue, C. Horizontal and vertical movements of humpback whales inform the use of critical pelagic habitats in the western South Pacific. Sci. Rep. 10, 1–13 (2020).Article 

Google Scholar 
Garrigue, C. et al. First assessment of interchange of humpback whales between Oceania and the East coast of Australia. J. Cetac. Res. Manage. 3, 269–274 (2011).
Google Scholar 
Steel, D. et al. Migratory connections between humpback whales from South Pacific breeding grounds and Antarctic feeding areas based on genotype matching. Int. Whal. Comm. (2008).Constantine, R., Russell, K., Gibbs, N., Childerhouse, S. & Baker, C. S. Photo-identification of humpback whales (Megaptera novaeangliae) in New Zealand waters and their migratory connections to breeding grounds of Oceania. Mar. Mam. Sci. 23, 715–720 (2007).Article 

Google Scholar 
Garland, E. C. et al. Humpback whale song on the southern ocean feeding grounds: implications for cultural transmission. PLoS ONE 8, e79422 (2013).ADS 
Article 

Google Scholar 
Garland, E. C. et al. Population structure of humpback whales in the western and central South Pacific Ocean as determined by vocal exchange among populations. Conserv. Biol. 29, 1198–1207 (2015).Article 

Google Scholar 
Cholewiak, D. M., Sousa-Lima, R. S. & Cerchio, S. Humpback whale song hierarchical structure: Historical context and discussion of current classification issues. Mar. Mam. Sci. 29, E312–E332. https://doi.org/10.1111/mms.12005 (2013).Article 

Google Scholar 
Payne, K., Tyack, P. & Payne, R. Progressive changes in the songs of humpback whales (Megaptera novaeangliae): a detailed analysis of two seasons in Hawaii. Communication and behavior of whales, 9–57 (1983).Payne, K. & Payne, R. Large scale changes over 19 years in songs of humpback whales in Bermuda. Ethology 68, 89–114 (1985).
Google Scholar 
Allen, J. A., Garland, E. C., Dunlop, R. A. & Noad, M. J. Cultural revolutions reduce complexity in the songs of humpback whales. Proc. R. Soc. B: Biol. Sci. 285, 20182088. https://doi.org/10.1098/rspb.2018.2088 (2018).Article 

Google Scholar 
Allen, J. A., Garland, E. C., Murray, A., Noad, M. J. & Dunlop, R. Using self-organizing maps to classify humpback whale song units and quantify their similarity. J. Acoust. Soc. Am. 142, 1943–1952 (2017).ADS 
Article 

Google Scholar 
Murray, A., Dunlop, R. A., Noad, M. J. & Goldizen, A. W. Stereotypic and complex phrase types provide structural evidence for a multi-message display in humpback whales (Megaptera novaeangliae). J Acoust Soc Am. 143, 980–994 (2018).ADS 
Article 

Google Scholar 
Garland, E. C. et al. Redefining western and central South Pacific humpback whale population structure based on vocal cultural exchange. (2013).Rekdahl, M. Humpback whale vocal communication: Use and stability of social calls and revolutions in the songs of east Australian whales. (2012).Templeton, C. N., Laland, K. N. & Boogert, N. J. Does song complexity correlate with problem-solving performance in flocks of zebra finches?. Anim. Behav. 92, 63–71 (2014).Article 

Google Scholar 
Boogert, N. J., Giraldeau, L.-A. & Lefebvre, L. Song complexity correlates with learning ability in zebra finch males. Anim. Behav. 76, 1735–1741 (2008).Article 

Google Scholar 
Winn, H. & Winn, L. The song of the humpback whale Megaptera novaeangliae in the West Indies. Mar. Biol. 47, 97–114 (1978).Article 

Google Scholar 
Girola, E., Noad, M. J., Dunlop, R. A. & Cato, D. H. Source levels of humpback whales decrease with frequency suggesting an air-filled resonator is used in sound production. The Journal of the Acoustical Society of America (In Review).Warren, V. E., Constantine, R., Noad, M., Garrigue, C. & Garland, E. C. Migratory insights from singing humpback whales recorded around central New Zealand. R. Soc. Open Sci. 7, 201084 (2020).ADS 
Article 

Google Scholar 
Garland, E. C. et al. Quantifying humpback whale song sequences to understand the dynamics of song exchange at the ocean basin scale. J. Acoust. Soc. Am. 133, 560–569. https://doi.org/10.1121/1.4770232 (2013).ADS 
Article 
PubMed 

Google Scholar 
Owen, C. et al. Migratory convergence facilitates cultural transmission of humpback whale song. R. Soc. Open Sci. 6, 190337 (2019).ADS 
Article 

Google Scholar 
Garland, E. C., Rendell, L., Lamoni, L., Poole, M. M. & Noad, M. Song hybridization events during revolutionary song change provide insights into cultural transmission in humpback whales. Proc. Natl. Acad. Sci. 114, 7822–7829 (2017).CAS 
Article 

Google Scholar 
Noad, M. & Cato, D. A combined acoustic and visual survey of humpback whales off southeast Queensland. Mem. Qld. Mus. 47, 507–523 (2001).
Google Scholar 
Spierings, M., de Weger, A. & Ten Cate, C. Pauses enhance chunk recognition in song element strings by zebra finches. Anim. Cogn. 18, 867–874 (2015).Article 

Google Scholar 
Doupe, A. J. & Kuhl, P. K. Birdsong and human speech: common themes and mechanisms. Annu. Rev. Neurosci. 22, 567–631 (1999).CAS 
Article 

Google Scholar 
Allen, J. A., Garland, E. C., Dunlop, R. A. & Noad, M. J. Network analysis reveals underlying syntactic features in a vocally learnt mammalian display, humpback whale song. Proc. R. Soc. B 286, 20192014 (2019).Article 

Google Scholar 
Barón Birchenall, L. Animal communication and human language: An overview. International Journal of Comparative Psychology 29 (2016).Noad, M. J. The use of song by humpback whales (Megaptera novaeangliae) during migration off the east coast of Australia (doctoral dissertation) Doctor of Philosophy thesis, University of Sydney, (2002).Catchpole, C. Song and female choice: good genes and big brains?. Trends Ecol. Evol. 11, 358–360. https://doi.org/10.1016/0169-5347(96)30042-6 (1996).CAS 
Article 
PubMed 

Google Scholar 
Nowicki, S., Hasselquist, D., Bensch, S. & Peters, S. Nestling growth and song repertoire size in great reed warblers: evidence for song learning as an indicator mechanism in mate choice. Proc. R. Soc. London B: Biol. Sci. 267, 2419–2424 (2000).CAS 
Article 

Google Scholar 
NOAA. Vol. 81 (ed National Marine Fisheries Service) 62260–62320 (Department of Commerce, Federal Register, 2016).Noad, M. J., Kniest, E. & Dunlop, R. A. Boom to bust? Implications for the continued rapid growth of the eastern Australian humpback whale population despite recovery. Popul. Ecol. 61(2), 198–209 (2019).Article 

Google Scholar 
Garrigue, C., Albertson, R. & Jackson, J. A. An anomalous increase in the New Caledonia humpback whales breeding sub-stock E2. Scientific Committee of the International Whaling Commission, Paper (2012).Garland, E. C. & McGregor, P. K. Cultural transmission, evolution, and revolution in vocal displays: Insights from bird and whale song. Front. Psychol. 11, 2387 (2020).Article 

Google Scholar 
Zandberg, L., Lachlan, R. F., Lamoni, L. & Garland, E. C. Global cultural evolutionary model of humpback whale song. Philos. Trans. R. Soc. B 376, 20200242 (2021).Article 

Google Scholar 
Crates, R. et al. Loss of vocal culture and fitness costs in a critically endangered songbird. Proc. R. Soc. B 288, 20210225 (2021).Article 

Google Scholar 
Garland, E. C., Garrigue, C. & Noad, M. J. When does cultural evolution become cumulative culture? A case study of humpback whale song. Philos. Trans. R. Soc. B 377, 20200313 (2022).Article 

Google Scholar 
Garland, E. C. et al. Improved versions of the Levenshtein distance method for comparing sequence information in animals’ vocalisations: tests using humpback whale song. Behaviour 149, 1413–1441 (2012).Article 

Google Scholar 
Garland, E. C. et al. The devil is in the detail: quantifying vocal variation in a complex, multileveled, and rapidly evolving display. J. Acoust. Soc. Am. 142, 460–472 (2017).ADS 
Article 

Google Scholar 
Rekdahl, M. L. et al. Culturally transmitted song exchange between humpback whales (Megaptera novaeangliae) in the southeast Atlantic and southwest Indian Ocean basins. R. Soc. Open Sci. 5, 172305 (2018).ADS 
Article 

Google Scholar 
Suzuki, R. & Shimodaira, H. Pvclust: an R package for assessing the uncertainty in hierarchical clustering. Bioinformatics 22, 1540–1542 (2006).CAS 
Article 

Google Scholar 
Sokal, R. R. & Rohlf, F. J. The comparison of dendrograms by objective methods. Taxon 11(2), 33–40 (1962).Article 

Google Scholar 
R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2015. URL http (R Core Development Team, 2016). More

These authors contributed equally: Gábor Pozsgai, Ibtissem Ben Fekih.State Key Laboratory of Ecological Pest Control for Fujian and Taiwan Crops, Institute of Applied Ecology, Fujian Agriculture and Forestry University, Fuzhou, 350002, ChinaGábor Pozsgai, Ibtissem Ben Fekih, Jie Zhang & Minsheng YouJoint international Research Laboratory of Ecological Pest Control, Ministry of Education, Fuzhou, 350002, ChinaGábor Pozsgai, Gábor L. Lövei & Minsheng YouCE3C – Centre for Ecology, Evolution and Environmental Changes, Azorean Biodiversity Group and Universidade dos Açores, Angra do Heroísmo, 9700-042, Azores, PortugalGábor PozsgaiInstitute of Environmental Microbiology, College of Resources and Environment, Fujian Agriculture and Forestry University, Fuzhou, 350002, ChinaIbtissem Ben Fekih & Christopher RensingBasic Forestry and Proteomics Research Center, College of Life Science, Fujian Provincial Key Laboratory of Haixia Applied Plant Systems Biology, Fujian Agriculture and Forestry University, Fuzhou, 350002, ChinaMarkus V. KohnenLaboratoire de Biologie et de Physiologie des Organismes, Faculté des Sciences Biologiques, Université des Sciences et de la Technologie Houari Boumediène, BP 32 El Alia, Alger, 16111, AlgeriaSaid AmraniDuna-Ipoly National Park Directorate, Költő u. 21, H-1121, Budapest, HungarySándor BércesJuhász-Nagy Pál Doctoral School, University of Debrecen, Egyetem tér 1, H-4032, Debrecen, HungarySándor BércesDepartment of Zoology, Plant Protection Institute, Centre for Agricultural Research, Nagykovácsi út 26-30, H-1029, Budapest, HungaryDávid FülöpFujian University Key Laboratory for Plant-Microbe Interaction, College of Plant Protection, Fujian Agriculture and Forestry University, Fuzhou, 350002, ChinaMohammed Y. M. JaberDepartment of Plant and Environmental Sciences, University of Copenhagen, Thorvaldsensvej 40, 1871, Frederiksberg C, DenmarkNicolai Vitt MeylingDepartment of Algology and Mycology Faculty of Biology and Environmental Protection, University of Łódź, Banacha 12/16, PL-90-237, Łódź, PolandMalgorzata Ruszkiewicz-MichalskaDepartment of Molecular Biotechnology and Microbiology, University of Debrecen, Egyetem tér 1, Debrecen, H-4032, HungaryWalter P. PflieglerFujian Provincial Key Laboratory of Insect Ecology, College of Plant Protection, Fujian Agriculture and Forestry University, Fuzhou, 350002, ChinaFrancisco Javier Sánchez-GarcíaÁrea de Biología Animal, Departamento de Zoología y Antropología Física, Facultad de Veterinaria, Universidad de Murcia, Murcia, 30100, SpainFrancisco Javier Sánchez-GarcíaDepartment of Agroecology, Aarhus University, Flakkebjerg Research Centre, Forsøgsvej 1, DK-4200, Slagelse, DenmarkGábor L. Lövei More

Effect of warming on physical and chemical conditionsThe water temperature in the warmed treatment was increased by 2.87 ± 0.20 °C in spring and 3.04 ± 0.08 °C in fall, compared to the control (Fig. 1a,b, Table 1). The average temperature in the control treatment, throughout the duration of the experiment, was about 4 °C cooler in spring (14.84 ± 0.03 °C) than in fall (19.01 ± 0.02 °C). In spring, the temperature naturally increased by approximately 4.19 °C from day (d) 10, until the end of the experiment, whereas it remained relatively constant in the fall experiment. It displayed higher diurnal variations in spring than in fall: over the course of the experiments, daily temperature variation ranged from 0.87 to 1.98 °C in spring and from 0.23 to 1.12 °C in fall. The average Daily Light Integral (DLI) in the control treatment, was almost twice as high in the spring experiment (7.93 ± 0.61 mol m−2 d−1) than during the fall (4.61 ± 0.52 mol m−2 d−1) (Fig. 1c,d). Warming did not significantly alter the DLI in fall (Table 1); however, the DLI could not be measured in the warm mesocosms in spring, owing to technical problems.Figure 1Time series of physical and chemical variables. Water temperature (a, b), Daily Light Integral (DLI, c, d), ammonium (NH4+, e, f), nitrates (NO3− + NO2−, g, h), orthophosphate (PO43−, i, j), silicate (SiO2, k, l) concentrations, and N/P ratio (m, n) over the course of the spring (a, c,e, g, i, k, m) and fall (b, d, f, h, j, l, n) experiments in the control (black) and the warmed (orange) treatments. Error bars represent range of observation for the two mesocosms per treatment in spring and the standard deviation for the three mesocosms per treatment in fall. Dotted lines represent the missing data on d10 of the fall experiment due to bad weather conditions. Due to technical difficulties, DLI could not be calculated in the warmed mesocosms of the spring experiment.Full size imageTable 1 Summary of the p values obtained by Repeated Measures Analyses Of VAriance (RM-ANOVA, with treatment as fixed factor and time as random factor) comparing physical parameters and nutrient concentrations in the warmed and control mesocosms.Full size tableNutrient concentrations were measured daily in all mesocosms (Fig. 1, Table 1). Ammonium (NH4+) concentrations were higher in spring than in fall in the controls (0.45 ± 0.08 µM, and 0.41 ± 0.05 µM, respectively). Ammonium concentrations were significantly different, between the control and warmed mesocosms, only at the end of the fall experiment (warmed with a mean of 0.58 ± 0.23 µM, between d11–17; and control with a mean of 0.21 ± 0.09 µM, between d11–17). In contrast, ammonium concentrations did not vary between the control and warmed mesocosms in spring (Table 1). Cohen’s effect size (d) was used to evaluate the magnitude of the effect of warming. Regarding ammonium concentrations, the values were ten times larger in spring than in fall.The nitrate + nitrite (NO3− + NO2−) concentrations in the control treatments were higher in the spring, compared to the fall experiment (0.71 ± 0.08 µM and 0.23 ± 0.02 µM, respectively). Moreover, warming had different effects, depending on the experiment. In spring, the nitrate + nitrite concentrations were significantly higher in the warmed mesocosms from d8 until the end of the experiment, with an average difference of 540.5% between the warmed and control mesocosms, corresponding to a very large d. In fall, nitrate + nitrite concentrations were significantly lower in the warmed mesocosms than in the control, with an average difference of 20.4%, and a medium-sized d.Similar to the nitrate + nitrite concentrations, the orthophosphate (PO43−) concentrations in the control treatment were higher in spring (0.55 ± 0.07 µM and 0.17 ± 0.01 µM, respectively) than in fall. The concentrations were negatively affected by warming throughout the spring experiment, with an average decrease of 9.3%. However, the largest negative effect of experimental warming was observed at the end of fall, with an average decrease of 16.7%, between d15 and d17.Contrary to the nitrate + nitrite and orthophosphate concentrations, the silicate (SiO2) concentrations in the control treatments were, on average, lower in spring (3.31 ± 0.18 µM and 10.42 ± 0.15 µM, respectively) than in fall. Warming had a significant positive effect during the second part of the spring experiment (from d10 to d17), with average concentrations being 27.8% higher in the warmed mesocosms, than in the control. In contrast, the strongest effect of warming on silicate concentrations was observed at the end of the fall experiment, when the silicate concentrations were significantly lower in the warmed mesocosms, by an average of 10.8%, between d15 and d17.The N/P ratio, calculated as the sum of nitrate, nitrite and ammonium concentrations divided by orthophosphate concentration, was 1.99 and 2.39 on average in the control treatment of the spring and fall experiments, respectively (Fig. 1m,n). It was significantly higher over the entire experiments in the warmed treatment by on average 191.5% and 58.8% in spring and fall, respectively (Table 1). In fall, the highest difference between treatments was seen during the second half of the experiment, when the ratio was significantly higher by 133.5% from day 11 to 17.Effects of warming on gross primary production and respiration rates derived from oxygen sensor dataIn the spring experiment, daily GPP varied between 0.19 ± 0.01 and 1.72 ± 0.15 gO2 m−3 d−1 in the control mesocosms (Fig. 2A). It increased during the first half of the experiment (d2–d10), then decreased toward the end of the experiment. In the fall experiment, the daily GPP was lower than what was observed in the control mesocosms in spring and varied between 0.12 ± 0.02 and 0.96 ± 0.16 gO2 m−3 d−1 (Fig. 2B). In spring, warming significantly reduced GPP by 50.9% over the entire experiment, while in fall, warming enhanced GPP by 21.1% over the entire experiment, 32.3% from d4 to d7, and 44.1% from d12 to d17 (Table 2). In spring, when GPP was normalized by the chl-a measured by the high-frequency sensors, it was not significantly different between the treatments, over the entire experiment (Fig. 2C, Table 2). However, it was significantly higher (138%) in the warmed treatment, during the second half of the experiment (d10– d17). In fall, GPP normalized by the chl-a was also significantly enhanced (12%) by warming (Fig. 2D, Table 2).Figure 2Plankton oxygen metabolism parameters. Gross Primary Production (GPP, A, B), GPP normalized by the chlorophyll-a fluorescence (C, D), Respiration (R, E, F), R normalized by the chlorophyll-a fluorescence (G, H), and GPP:R ratio (I, J) in the control (black) and warmed (orange) treatments. Error bars represent range of observation for the two mesocosms per treatment in spring and the standard deviation for the three mesocosms per treatment in fall. In the spring experiment, GPP: Chl-a and R: Chl-a could not be estimated on d1 and d2.Full size imageTable 2 Summary table of the p values and the F-values obtained with the RM-ANOVA (with treatment as fixed factor and time as random factor) comparing the chl-a fluorescence, µ, l, the µ:l ratio, and pigment concentrations in the warmed and in the control treatments over the entire spring and fall experiments or over specific periods defined after trends observed in the data.Full size tableDaily R varied between 0.27 ± 0.02 and 1.92 ± 0.20 gO2 m−3 d−1 in the spring control mesocosms (Fig. 2E). Similar to the daily GPP, it increased during the first half of the experiment (d2– d10), with a strong increase between d8 and d10, before decreasing slowly until the end of the experiment. In the fall experiment, the daily R was lower than in spring, varying from 0.19 ± 0.02 and 1.09 ± 0.14 gO2 m−3 d−1, in the control mesocosms (Fig. 2F). Warming significantly reduced the daily R by an average of 47.9% in spring, while no significant differences were found in fall (Table 2). During both experiments, when daily R was normalized by chl-a, it was not significantly different between treatments over the entire experimental period (Figs. 2G,H), but it was significantly enhanced by warming during the second half of the experiments, by 172% and 49.6%, from d10–17 in spring, and d11–17 in fall, respectively (Table 2).The GPP:R ratio was on average 1.01 and 1.08 in the spring and fall control treatments, respectively (Fig. 2I,J). Consequently, because warming decreased GPP and R to a similar extent in spring, it did not significantly change the GPP:R ratio. Warming significantly increased GPP:R, by an average of 32% in fall (Table 2).Effects of warming on phytoplankton biomass (chlorophyll-a), growth, and loss rates derived from the chlorophyll-a sensor dataThe chl-a fluorescence data was measured using high-frequency sensors, which were inter-calibrated before and after the experiments, and were corrected by the chl-a concentration measured daily by HPLC (see “Methods”). It is hereafter referred to as chl-a. In the spring experiment, the daily chl-a was 5.28 ± 0.21 µg L−1 in the control mesocosms (Fig. 3a,c). A phytoplankton bloom dynamic was observed, with increasing concentrations from d2 to d10, reaching a maximum value of 8.62 ± 0.15 µg L−1, and decreasing concentrations from d10 to d17. The average daily chl-a was lower in fall than in the spring experiment (4.30 ± 0.59 µg L−1) (Fig. 3b,d), and displayed a relatively flat dynamic during the entire experiment, with maximum values on d8 (5.53 ± 0.58 µg L−1).Figure 3Phytoplankton chlorophyll-a, growth and loss rates. High-frequency chlorophyll-a data, uncorrected for Non Photochemical Quenching (NPQ) (a, b), daily average chlorophyll-a data corrected for the NPQ(c, d), phytoplankton growth rate (µ, e, f), loss rate (l, g, h), and µ:l ratio (i, j) in the control (black) and warmed (orange) treatments. Error bars represent range of observation for the two mesocosms per treatment in spring and the standard deviation for the three mesocosms per treatment in fall. In the spring experiment, µ and l could not be estimated on d1 and d12 and, for the latter, the missing data are represented as dotted lines.Full size imageWarming significantly reduced chl-a in both experiments (Table 2): an average of 69.5% from d5 to the end of the spring, and 31.7% from d8 to 15, in the fall experiment. Conversely, warming significantly enhanced chl-a concentrations at the beginning of the fall experiment (19.4% between d2 and d6). Generally, the magnitude of the effect was larger in spring than in fall (Table 2).In the control treatment, µ was higher in spring than in fall (0.44 ± 0.04 d−1 and 0.32 ± 0.05 d−1, respectively; Fig. 3e,f). During both seasons, the maximum µ was observed during the first half of the experiment (spring d7, 0.99 ± 0.01 d−1; fall d4, 0.61 ± 0.03 d−1). Warming enhanced µ by an average of 18.3% and 28.1%, over the entire spring and fall experiments, respectively, and by an average of 56.8% and 50.9%, respectively, from d8 until the end of the experiment (Table 2). The effect size was higher in fall than in spring (Table 2). However, contrary to the general trend of the entire experiment, during spring, warming significantly reduced µ during the first part of the experiment (d2–d7), with an 18.8% mean difference between the treatments.In contrast to µ, l was almost similar between the seasons, with average values of 0.39 ± 0.04 d−1 and 0.40 ± 0.07 d−1, in the control treatments for spring and fall, respectively (Fig. 3g,h). In the spring experiment, warming had a positive effect on the mean l across the study period (37.1%), and even more from d8 to d17 (59.1%), which was larger than the positive effect found for µ. The effect size of warming was not as large in fall, and l was significantly higher in the warmed treatment, although only in the middle of the experiment (20.4% from d7 to d11).When comparing µ and l, the results showed that in spring, µ was higher than l in the control treatment, during the first part of the experiment (d2–d9), and lower during the latter half of the experiment (d10–d17). Warming significantly decreased the µ:l ratio by 28.9%, during the first half of the experiment (D 3–8, Fig. 3i, Table 2), whereas no significant effect was observed in the rest of the experiment. In the fall control treatment, the µ:l ratio was generally lower than that of the spring control (Fig. 3j). Contrary to what was observed in the spring warming, this ratio significantly increased by an average of 92.9%, in the second half of the experiment (d11–d17, Table 2).Effects of warming on phytoplankton pigment concentrationsPhytoplankton pigment composition varied between seasons (Fig. 4). In the spring control treatment, the predominant pigments were fucoxanthin and 19′-hexanoyloxyfucoxanthin (19′-HF), which are mostly associated with diatoms (1.14 ± 0.10 µg L−1) and prymnesiophytes (19′-HF, 2.91 ± 0.14 µg L−1)23,24, respectively (Figs. 4A,B,G,H). The other pigments that were present included peridinin (0.18 ± 0.01 µg L−1), Chl-b (0.14 ± 0.01 µg L−1), zeaxanthin (0.08 ± 0.01 µg L−1), and the specific accessory pigment prasinoxanthin (0.06 ± 0.01 µg L−1), which are associated with dinoflagellates, green algae, cyanobacteria, and prasinophytes, respectively (Figs. 4C–F,I,J)23,24,25.Figure 4Phytoplankton pigment concentrations. Daily pigment concentrations (µg L−1) in the control (black) and warmed (orange) treatments for the spring (A–F) and fall (G–J) experiments. Error bars represent range of observation for the two mesocosms per treatment in spring and the standard deviation for the three mesocosms per treatment in fall. Dotted lines represent the missing data on d10 of the fall experiment due to bad weather conditions. (A, G) fucoxanthin; (B, H) 19′-hexanoyloxyfucoxanthin; (C, I) zeaxanthin; (D, J) chlorophyll-b; (E) peridinin, and (F) prasinoxanthin. Corresponding phytoplankton functional groups are indicated in parentheses.Full size imageIn the fall control treatment, the dominant pigments were the cyanobacteria-associated zeaxanthin (1.78 ± 0.22 µg L−1), the diatom-associated fucoxanthin (1.07 ± 0.27 µg L−1), the green algae-associated Chl-b (0.69 ± 0.14 µg L−1), and the prymnesiophyte-associated 19′-HF (0.68 ± 0.16 µg L−1). Among the main pigments that were identified in the spring experiment, peridinin and prasinoxanthin were either not detected or detected at negligible concentrations in the fall experiment, whereas lutein was detected in fall but not in spring (data not shown).Warming had seasonal effects on pigment concentrations (Table 2). In the spring experiment, warming had a large and significant negative effect on 19′-HF and zeaxanthin concentrations, with mean concentrations decreasing by 75.4% and 75.2%, respectively. Conversely, warming had moderately significant positive effects on peridinin concentration, which increased by an average of 101%.In the fall experiment, warming had a significant negative effect on Chl-b concentration, which decreased by 19.5%, and on zeaxanthin concentration, which significantly decreased in the middle of the experiment (43.4% from d11 to d15). In contrast, a significant positive effect was observed on fucoxanthin concentration, which increased by 210.7%, during the second part of the experiment (between d13 and d17).Relationships between plankton processes, pigment concentrations and environmental parametersPrincipal component analyses (PCA) were used to project plankton processes, pigment concentrations and environmental parameters in a multidimensional space in order to illustrate relationships among variables in both experiments (Fig. 5). For both experiments, GPP and R were clustered together along the first PCA axis, although they appeared closer in spring than in fall (Fig. 5A,B). Conversely, µ was close to ammonium for both experiments, to silicate in spring and to nitrate and nitrite in fall; and l was part of this cluster in spring but not in fall. Concerning phytoplankton pigment composition, in spring, zeaxanthin, associated with cyanobacteria, and 19′-HF, associated with prymnesiophytes, were part of a group together with temperature (Fig. 5C). Similarly, prasinoxanthin, associated with prasinophytes, and Chl-b, associated with green algae, were grouped with DLI and orthophosphate. Finally, peridinin, associated with dinoflagellates, and silicate were clustered together and opposed to fucoxanthin, which is associated with diatoms. In fall, zeaxanthin and Chl-b, representing cyanobacteria and green algae, were part of a group opposed to N-nutrients and temperature, while fucoxanthin was opposed to DLI, silicate and orthophosphate (Fig. 5D).Figure 5Principal component analyses (PCA) of logarithm response ratio (LRR) of plankton processes (A, B) and pigment concentrations (C, D) with environmental parameters for the spring (A, C) and fall (B, D) experiments. GPP: Gross Primary Production, R: Respiration, µ: Phytoplankton growth rate, l: Phytoplankton loss rate, Chl-b: Chlorophyll-b, 19′-HF: 19′-Hexanoyloxyfucoxanthin, Fuco: Fucoxanthin, Prasino: Prasinoxanthin, Zea: Zeaxanthin, DLI: Daily Light Integral.Full size imageTo evaluate specific relationships between phytoplankton processes, environmental variables, and phytoplankton community composition, ordinary least squares linear relationships were assessed for the effects of warming (expressed as the logarithmic response ratio) on GPP, R, µ, and l, nutrient concentrations, DLI, and pigment concentrations (Fig. 6). A significant positive relationship was found between the effects of warming on GPP versus R, and µ versus l (Fig. 6A,B). Moreover, the effect of warming on µ was positively and linearly related to the effects of warming on ammonium in both seasons (Fig. 6C), and to nitrate + nitrite concentrations in spring (Fig. 6D). There was no relationship between the effects of µ on pigment concentrations in the spring experiment, but its effects were positively correlated with the diatom-associated pigment fucoxanthin in fall (Fig. 6E). Similarly, R was positively correlated with fucoxanthin in fall (Fig. 6F). In contrast, significant negative relationships were found between the effects of warming on µ, Chl-b, and zeaxanthin, the pigments associated with green algae and cyanobacteria, respectively (Fig. 6G,H). Similarly, the effect of warming on R was negatively correlated to orthophosphate (Fig. 6I), and the effect on GPP to nitrate + nitrite and peridinin concentrations (Fig. 6J,K).Figure 6Linear relationships between the effect of warming on plankton processes, environment variables and pigment concentrations. Ordinary least squares linear relationships between the effect of warming, expressed as the log response ratio, on GPP, R, µ, and l, and the effect of warming on environmental and pigment variables for the spring (blue circles) and fall (green squares) experiments. Relationships were individually assessed for each experiment. Only statistically significant relationships (p  More

Choice of biowaste treatment plants and sample identifiersCompost samples were collected from four central municipal biowaste treatment plants (denominated as #1 to #4) in Baden-Wurttemberg, Germany (Table 1). All plants used a state-of-the-art two-stage biowaste treatment process comprising of (a) anaerobic digestion/biogas production and (b) subsequent composting of the solid digestate to produce a high-quality mature compost sold for direct use as fertilizer in agriculture. The composts were regularly analyzed by an independent laboratory for quality and residual contamination and consistently fulfilled the quality requirements of the label RAL-GZ 251 Gütezeichen Kompost of the German Bundesgütegemeinschaft Kompost e.V. (www.gz-kompost.de). Plants #1 and #3 produce in addition a liquid fertilizer, which is separated from the solid digestate at the end of stage a) by press filtration and which is also intended for direct use on agricultural soil (replacement of liquid manure). In case of plants #1, #3, and #4 up to 25 wt% of shrub/tree cuttings were added to the solid digestate for composting. All plants used sieving (typically with a 12 or a 20 mm mesh) at the end of the process to assure the necessary purity of their finished composts. Whenever technically possible, we as well took samples of the pre-compost immediately before this final sieving step to evaluate its contribution to the removal of residual BPD fragments. For analysis, composts were passed consecutively through two sieves with mesh sizes of 5 mm and 1 mm, yielding two fragment preparations for IR-analysis namely a > 5 mm fraction corresponding to the contamination by residual “macroplastic” (5 mm is a commonly used upper size limit for “microplastic”, anything larger is macroplastic) and a 1–5 mm fraction corresponding to the regulatory relevant residual contamination by microplastic. The lower limit of 1 mm rather than 2 mm was chosen in anticipation of the expected changes in regulation, where the replacement of the 2 mm limit by a 1 mm limit is imminent.Table 1 Technical data of the investigated plants and incidence of BDP fragments in the sampled composts.Full size tableOccurrence of plastic fragments  > 1 mm in the sampled compostsComposting times of 5–9 weeks were used in the investigated plants (Table 1), which is shorter than the 12 weeks indicated in EN 13432 for the 90% disintegration of a compostable plastic material, but a realistic time span for state-of-the-art technical waste treatment. Since we were not in a position to estimate the quantity of BDP entering the plants, since for technical reasons we were unable to obtain a representative sample, we cannot say, whether any residual BDP detected by us in the finished composts was due to a yet incomplete disintegration process or whether it corresponds to the 10% material still permissible by EN 13432 even after the full composting step. However, in 7 out of the 12 sampled composts and pre-composts fragments with chemical signatures corresponding to the BDPs poly (lactic acid) (PLA) and poly (butylene-adipate-co-terephthalate) (PBAT) were identified in the > 5 mm and/or the 1–5 mm sieving fractions using FTIR analysis3 (Fig. 1; Table 1). All recovered fragments appeared to stem from foils, bags or packaging, since they were thin compared to their length and width (see Suppl Figure S1 for typical examples). Fragments with overlapping signatures, most likely PBAT/PLA mixtures or blends, were also found (see Suppl Figure S2 for the interpretation of the spectra). In addition, the recorded BDP fragment spectra (Fig. 1A) showed high similarity to the FTIR spectra of commercial compostable bags sold in the vicinity of the biowaste treatment plants (Fig. 1B), which together with the geometry of the recovered fragments led us to assuming that the majority of the BDP entered the biowaste in the form of such bags.Figure 1FTIR spectra of BDP fragments from composts and commercial bags. (A) BDP fragments recovered from the composts and (B) the commercial compostable bags. Fragments were coded as follows: p or f for pre-compost or finished compost, followed by the plant number (#1 to #4), an indication of the size fraction ( > 5 mm or 1–5 mm) in which the fragment was found, and finally, the fragment number. Fragment F#1_5mm_4 therefore represents the 4th fragment collected in the  > 5 mm size fraction from the finished compost of plant number 1. Bags were arbitrarily numbered 1–10, see Suppl Table S1 for supplier information. The spectra (in grey) of the reference materials for PLA and PBAT are given as basis for the interpretation. Spectra in red refer to test samples consisting only of PBAT, while those in blue indicate samples composed of PBAT/PLA mixtures.Full size imageThe BDP fragments were found alongside fragments of commodity plastics (mostly PE) in all cases. Finished composts tended to contain fewer and smaller fragments than the corresponding pre-composts. The final sieving of the pre-composts to prepare the finished composts hence appears to be quite effective in removing such fragments, in particular those from the > 5 mm size fraction (Table 1) and for that reason has become state-of-the-art in preparing quality composts (contamination by plastic fragments > 2 mm of less than 0.1 wt%). Given that the size of the fragments is a crucial factor regarding ecological risk, we analyzed the sizes (length Î width) of the BDP fragments in comparison to that of the plastic fragments with signatures of commodity plastics such as PE (Fig. 2). BDP fragments found in a given compost sample tended to be smaller than the fragments stemming from non-BDP materials, which may indicate that BDPs degrade faster or tend to disintegrate into tinier particles than commodity plastics. This may also explain why in the compost from plant #2, no BDP fragments were found in the particle fraction retained by the 5 mm sieve ( > 5 mm fraction), while 19 such particles were found in the fraction then retained by the 1 mm sieve (1–5 mm fraction). Interestingly, plant #2 is the only one included in our study that uses no mechanical breakdown of the incoming biowaste. This reduces the mechanical stress on the incoming material. Mechanical stress can alter the properties of plastic foils such as the crystallinity whereby crystallinity has been shown to influence the biological degradation of BDP such as PLA7.Figure 2Size distribution of plastic fragments  > 1 mm. (A) Fragments found in the finished compost from plant #1, (B) in the finished compost from plant #2, and (C) in the pre-compost from plant #3. For reasons of statistical relevance, only samples containing more than 20 BDP fragments per kg of compost were included in the analysis.Full size imageMaterial characteristics of BDP fragments in comparison to those of commercial biodegradable bagsIn order to verify whether the BDP fragments recovered from the composts differed from the compostable bags in any parameter with possible relevance for biodegradation and environmental impact16, the physico-chemical properties of bags and fragments were studied in detail. Since we wanted to have a maximum of information of the BDP fragments, size/weight was a limiting factor in selecting fragments for analysis. Fragments of at least 1 mg were required for the FT-IR analysis. 5 mg-fragments could be analyzed in addition by 1H-NMR, while the full set of analytics (FT-IR, 1H-NMR, and DSC) required at least 10 mg of sample.For insight into the chemical composition, 1H-NMR spectra of the commercial bags and all suitable BDP fragments were compared (Fig. 3). In case of material mixtures and blends, the 1H-NMR analysis allows quantification of the PBAT/PLA weight ratio in the materials and also of the ratio of the butylene terephthalate (BT) and butylene adipate (BA) units in the involved PBAT polyesters.Figure 31H NMR spectra of BDP fragments from composts and commercial bags. (A) BDP fragments recovered from the composts and (B) the commercial compostable bags. Fragments were coded as follows: p or f for pre-compost or finished compost, followed by the plant number (#1 to #4), an indication of the size fraction ( > 5 mm or 1–5 mm) in which the fragment was found, and finally, the fragment number. Bags were arbitrarily numbered 1–10, see Suppl Table S1 for supplier information. The spectra (in grey) of the reference materials for PLA and PBAT are given as basis for the interpretation. Spectra in red refer to test samples consisting only of PBAT, while those in blue indicate samples composed of PBAT/PLA mixtures. (C) Chemical structures of PLA and PBAT, chemical shifts of the protons are assigned as indicated in the reference spectra in (B).Full size imageThe 1H-NMR spectra corroborate the FTIR measurements in that all investigated commercial bags were made from PBAT/PLA mixtures of varied composition (Table 2). By comparison, some of the fragments, for instance, f#1_5mm_4, appeared to consist of only PBAT. Other fragments, e.g., f#1_1mm_9, were mixtures of PLA and PBAT (Table 2). However, even in the case of PBAT/PLA mixtures, the average PBAT content tended to be higher in the fragments than in the bags, while the BT/BA monomer ratio in the respective PBATs, was also significantly higher in the fragments than in the bags. If we assume the fragments to stem from similar compostable bags as the ones included in our comparison, this would mean that during composting of such a bag, the PLA degrades more quickly than the PBAT, whereas within a given PBAT polyester, the BA unit is more easily degraded than the BT unit. Evidence can indeed be found in the pertinent literature that PLA has faster biodegradation kinetics than PBAT, while BT is more resistant to mineralization than BA17,18.Table 2 Composition of commercial compostable bags and BDP fragments recovered from the composts as analyzed by 1H-NMR.Full size tableNext, differential scanning calorimetry (DSC) was used to analyze fragments compared to commercial bags in regard to the presence of amorphous vs. crystalline domains, a parameter expected to affect biodegradation kinetics and therefore the putative environmental impact of the produced microplastic16 upon release into the environment with the composts. Whereas amorphous domains show glass transition, crystalline domains show melting, both of which can be discerned by the respective phase transition enthalpy in the DSC curves (Fig. 4).Figure 4DSC curves of BDP fragments and compostable bags #1 and #7. Curves for the reference materials (in grey) for PLA and PBAT are given for comparison. Curves were recorded during the first heating run (temperature range: − 50 °C to 200 °C, heating rate: 10 °C min−1). (A) and (B) curves in red refer to test samples consisting only of PBAT, while those in blue indicate samples composed of PBAT/PLA mixtures. Fragments were coded as follows: p or f for pre-compost or finished compost, followed by the plant number (#1 to #4), an indication of the size fraction ( > 5 mm or 1–5 mm) in which the fragment was found, and finally, the fragment number.Full size imageThe curve for the reference PBAT shows a glass transition temperature (Tg) of − 29 °C and a broad melting range between 100 and 140 °C for the crystalline domains, while that of the PLA reference shows a glass transition temperature of 58 °C and a narrower melting peak between 144 °C and 162 °C. The curve for commercial bag #1, which had a comparatively high PLA content, shows a pronounced melting peak in the expected range; the same is the case for fragment p#3_5mm_1 and to a lesser extent for fragment p#3_5mm_9, two fragments, which also have high PLA contents. The DSC curves of the other fragments and bag #1 are undefined in comparison, which is due to their high PBAT content. According to the DSC curves, most of the investigated materials are semicrystalline, i.e., contain both amorphous (glass transition) and crystalline (melting) domains. However, the DCS data alone allow only a qualitative discussion of the differences between fragments and bags.To obtain quantitative data on the crystallinity differences, wide angle X-ray scattering (WAXS) spectra were recorded. WAXS requires fragments at least 3 cm long, which restricted the number of fragment samples to three, all of which were found in pre-compost samples. The corresponding curves are shown in Fig. 5A–C. The spectra of the commercial biodegradable bags are shown in Suppl Figure S3. Foils were in addition prepared by heat pressing from the reference materials for PLA and PBAT in order to include them into the WAXS measurements (Fig. 5D). While the foils produced from the PBAT reference material produced crystallinity peaks at 16.2°, 17.3°, 20.4°, 23.2°, and 24.8°, the foil prepared from the PLA reference material showed only an amorphous halo at 15.5° and 31.5°, which is in accordance with values published in the literature19. A more pronounced crystallinity peak was obtained in the case of an additionally annealed PLA foil.Figure 5WAXS curves with Lorenz fitting for (A) fragment p#3_5mm_1, (B) fragment p#3_5mm_9, and (C) fragment p#4_5mm_2. (D) WAXS curves for foils produced from the PBAT and PLA reference materials; the percent values indicate the crystallinity. The dash lines are the fitting peak curves for the XRD spectrum. Crystallinity can be obtained by dividing the integration area of the fitted peaks by the integration area of the entire spectrum. Fragments were coded as follows: p or f for pre-compost or finished compost, followed by the plant number (#1 to #4), an indication of the size fraction ( > 5 mm or 1–5 mm) in which the fragment was found, and finally, the fragment number.Full size imageIn case of the fragments and bags, the peaks of PLA and PBAT overlapped to some extent in the WAXS spectra, but by conducting Lorenz fitting using Origin software, the overall crystallinity could be calculated as follows:$$chi = { 1}00% , *{text{ Aa}}/left( {{text{Aa }} + {text{ Ac}}} right)$$where χ is the crystallinity and Aa and Ac represent the areas of the amorphous and crystalline peaks.Using this equation, crystallinities of 55% (fragments p#3_5mm_1), 34% (p#3_5mm_9), and 34% (p#4_5mm_2) were calculated for the fragments. The foils prepared in house for the reference materials had similar crystallinities (43% in case of the annealed PLA foil and 26% of the PBAT foil), while the simple PLA foil was amorphous. By comparison, for eight of the commercial bags, crystallinities in the range from 1% to 7% were calculated, whereas these values were 14% and 15% for the remaining two bag types (Suppl Figure S3).The high crystallinity of the larger fragments recovered from the pre-compost samples suggests that crystalline domains of BDP materials may indeed disintegrate more slowly than the amorphous ones, as prior studies on microbial biodegradation have suggested7,8. Admittedly, such large fragments per se would not enter the environment, since the final sieving step used to prepare the finished composts is quite efficient at removing them. However, it is tempting to extrapolate that residual BDP in general are remnants of the more crystal domains of the original material, even though experimental proof of this assumption is at present not possible. 10 wt% of a BDP bag is allowed to remain after standard composting. It is usually assumed that any such residues continue to degrade with comparable speed. However, should these residues correspond to the more crystalline domains, rather than degrading with similar speed as the bulk material, the more crystalline fragments can be expected to persist for a much longer and at present unpredictable length of time in the environment, e.g. when applied to the soil with the composts; in particular, when they are also enriched in PBAT and BT units as suggested by our analysis of the chemical composition. Data from the use of biodegradable foils in agriculture show that the degradation in the environment may take years20. Altogether this may have unforeseen economic and environmental consequences, especially when considering the high fraction of BDP fragments < 5 mm. Putative consequences include changes in soil properties, the soil microbiome and therefore in plant performance21, a factor indispensable for worldwide nutrition.Residues of BDP fragments  1 mm were found in the collected LF samples. This is hardly surprising, given that the LF is produced by press filtration of the digestate after the anaerobic stage. Such a filtration step can be expected to retain fragments > 1 mm in the produced filter cake, which goes into the composting step, leaving the filtrate, i.e. the LF, essentially free of such particles. Anaerobic digestion is currently not assumed to contribute significantly to the degradation of BDP17,22, but the process conditions (mixing, pumping) may promote breakdown of larger fragments, particularly when additives such as plasticizers23 leach out of the material.Since the residual solids content of the LF is low (plant #1: 8.6 wt%, plant #3: 5.8 wt%), a combination of enzymatic-oxidative treatment and µFTIR imaging originally developed for environmental samples from aqueous systems24,25 could be adapted for the analysis (size and chemical signature) of particles in the LF down to a size of 10 µm. The corresponding data are compiled in Table 3. In all cases, residual fragments from PBAT-based polymers represented the dominant plastic fraction in the investigated samples; i.e. approximately 53% of all plastic particles in the LF from plant #1 (11,520 BDP particles per liter) and 65% in the case of plant #3 (12,480 BDP particles per liter). Liquid manure is applied several times a year to fields at a concentration of 2–3 L m−2. According to our analysis > 20,000 BDP microparticles of a size ranging from 10 µm to 500 µm enter each m2 of agricultural soil whenever LF is applied on agricultural surfaces.Table 3 Microplastic fragments (BDP/all) found per liter of liquid fertilizer.Full size tableDue to the complexity of the matrix, a similar analysis of individual plastic fragments  1 mm. Six compost samples representing the more contaminated ones based on the content of fragments > 1 mm, namely, f#1, f#2, p#3, f#3, p#4 and f#4 (nomenclature: f or p for finished or pre-compost, followed by plant number), were extracted with a 90/10 vol% chloroform/methanol mixture. The amounts of PBAT and PLA in the obtained extracts were then quantified via 1H-NMR (Table 4). Briefly, the intensity of characteristic signals in the extract spectra of the compost samples (see Suppl Figure S4) were compared to peak intensities produced by calibration standards of the pure polymer dissolved at a known concentration in the chloroform/methanol. All samples and standards were normalized using the 1,2-dichloroethan signal at 3.73 ppm as internal standard. See also Suppl Figure S5 for an exemplification of the quantification of the PBAT/PLA ratios. Based on the amounts of PBAT and PLA extracted from a known amount of compost, the total mass concentration (wt% dry weight) of these polymers in the composts was calculated.Table 4 Evidence of PBAT and PLA residues caused by fragments  2 mm. Moreover, residues of PBAT and PLA were found in all investigated compost samples, including the finished compost from plant #4, which had shown no contamination by larger BPD fragments (Table 1). The pre-compost from that plant had shown a few contaminating BDP fragments in the > 5 mm fraction. However, in regard to the fragments More