Philippa Kaur

Agostini, S. et al. Ocean acidification drives community shifts towards simplified non-calcified habitats in a subtropical–temperate transition zone. Sci. Rep. 8, 11354 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Walther, G.-R. et al. Ecological responses to recent climate change. Nature 416, 389–395 (2002).Article 
CAS 
PubMed 

Google Scholar 
Gilbert, B. & Levine, J. M. Ecological drift and the distribution of species diversity. Proc. Biol. Sci. 284, 20170507 (2017).PubMed 
PubMed Central 

Google Scholar 
White, E. P. et al. A comparison of the species–time relationship across ecosystems and taxonomic groups. Oikos 112, 185–195 (2006).Article 

Google Scholar 
Sax, D. F. & Gaines, S. D. Species invasions and extinction: the future of native biodiversity on islands. Proc. Natl Acad. Sci. USA 105, 11490–11497 (2008).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Thompson, J. N. Rapid evolution as an ecological process. Trends Ecol. Evol. 13, 329–332 (1998).Article 
CAS 
PubMed 

Google Scholar 
Post, D. M. & Palkovacs, E. P. Eco-evolutionary feedbacks in community and ecosystem ecology: interactions between the ecological theatre and the evolutionary play. Philos. Trans. R. Soc. Lond. B 364, 1629–1640 (2009).Article 

Google Scholar 
Reznick, D. N. & Travis, J. Experimental studies of evolution and eco-evo dynamics in guppies (Poecilia reticulata). Annu. Rev. Ecol. Evol. Syst. https://doi.org/10.1146/annurev-ecolsys-110218-024926 (2019).Hendry, A. P. Eco-evolutionary Dynamics (Princeton Univ. Press, 2020).Schoener, T. W. The newest synthesis: understanding the interplay of evolutionary and ecological dynamics. Science 331, 426–429 (2011).Article 
CAS 
PubMed 

Google Scholar 
Yoshida, T., Jones, L. E., Ellner, S. P., Fussmann, G. F. & Hairston, N. G. Rapid evolution drives ecological dynamics in a predator–prey system. Nature 424, 303–306 (2003).Article 
CAS 
PubMed 

Google Scholar 
Hansen, S. K., Rainey, P. B., Haagensen, J. A. J. & Molin, S. Evolution of species interactions in a biofilm community. Nature 445, 533–536 (2007).Article 
CAS 
PubMed 

Google Scholar 
Hillesland, K. L. & Stahl, D. A. Rapid evolution of stability and productivity at the origin of a microbial mutualism. Proc. Natl Acad. Sci. USA 107, 2124–2129 (2010).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Turcotte, M. M., Reznick, D. N. & Hare, J. D. The impact of rapid evolution on population dynamics in the wild: experimental test of eco-evolutionary dynamics. Ecol. Lett. 14, 1084–1092 (2011).Article 
PubMed 

Google Scholar 
Lawrence, D. et al. Species interactions alter evolutionary responses to a novel environment. PLoS Biol. 10, e1001330 (2012).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Celiker, H. & Gore, J. Clustering in community structure across replicate ecosystems following a long-term bacterial evolution experiment. Nat. Commun. 5, 4643 (2014).Article 
CAS 
PubMed 

Google Scholar 
Andrade-Domínguez, A. et al. Eco-evolutionary feedbacks drive species interactions. ISME J. 8, 1041–1054 (2014).Article 
PubMed 

Google Scholar 
Reznick, D. Hard and soft selection revisited: how evolution by natural selection works in the real world. J. Hered. 107, 3–14 (2016).Article 
PubMed 

Google Scholar 
Matthews, B., Aebischer, T., Sullam, K. E., Lundsgaard-Hansen, B. & Seehausen, O. Experimental evidence of an eco-evolutionary feedback during adaptive divergence. Curr. Biol. 26, 483–489 (2016).Article 
CAS 
PubMed 

Google Scholar 
Harcombe, W. R., Chacón, J. M., Adamowicz, E. M., Chubiz, L. M. & Marx, C. J. Evolution of bidirectional costly mutualism from byproduct consumption. Proc. Natl Acad. Sci. USA 115, 12000–12004 (2018).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Preussger, D., Giri, S., Muhsal, L. K., Oña, L. & Kost, C. Reciprocal fitness feedbacks promote the evolution of mutualistic cooperation. Curr. Biol. https://doi.org/10.1016/j.cub.2020.06.100 (2020).Adamowicz, E. M., Muza, M., Chacón, J. M. & Harcombe, W. R. Cross-feeding modulates the rate and mechanism of antibiotic resistance evolution in a model microbial community of Escherichia coli and Salmonella enterica. PLoS Pathog. 16, e1008700 (2020).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Rodríguez-Verdugo, A. & Ackermann, M. Rapid evolution destabilizes species interactions in a fluctuating environment. ISME J. 15, 450–460 (2021).Article 
PubMed 

Google Scholar 
Barber, J. N. et al. The evolution of coexistence from competition in experimental co-cultures of Escherichia coli and Saccharomyces cerevisiae. ISME J. 15, 746–761 (2021).Article 
CAS 
PubMed 

Google Scholar 
Hart, S. F. M., Chen, C.-C. & Shou, W. Pleiotropic mutations can rapidly evolve to directly benefit self and cooperative partner despite unfavorable conditions. eLife 10, e57838 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Kokko, H. et al. Can evolution supply what ecology demands? Trends Ecol. Evol. 32, 187–197 (2017).Article 
PubMed 

Google Scholar 
Nuismer, S. Introduction to Coevolutionary Theory (Macmillan Learning, 2017).Stoltzfus, A. & McCandlish, D. M. Mutational biases influence parallel adaptation. Mol. Biol. Evol. 34, 2163–2172 (2017).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Payne, J. L. et al. Transition bias influences the evolution of antibiotic resistance in Mycobacterium tuberculosis. PLoS Biol. 17, e3000265 (2019).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Storz, J. F. et al. The role of mutation bias in adaptive molecular evolution: insights from convergent changes in protein function. Philos. Trans. R. Soc. Lond. B 374, 20180238 (2019).Article 
CAS 

Google Scholar 
Gomez, K., Bertram, J. & Masel, J. Mutation bias can shape adaptation in large asexual populations experiencing clonal interference. Proc. Biol. Sci. 287, 20201503 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Venkataram, S., Monasky, R., Sikaroodi, S. H., Kryazhimskiy, S. & Kacar, B. Evolutionary stalling and a limit on the power of natural selection to improve a cellular module. Proc. Natl Acad. Sci. USA 117, 18582–18590 (2020).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Hom, E. F. Y. & Murray, A. W. Niche engineering demonstrates a latent capacity for fungal–algal mutualism. Science 345, 94–98 (2014).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Wolfe, B. E. & Dutton, R. J. Fermented foods as experimentally tractable microbial ecosystems. Cell 161, 49–55 (2015).Article 
CAS 
PubMed 

Google Scholar 
Chacón, J. M., Hammarlund, S. P., Martinson, J. N. V., Smith, L. B. & Harcombe, W. R. The ecology and evolution of model microbial mutualisms. Annu. Rev. Ecol. Evol. Syst. 52, 363–384 (2021).Article 

Google Scholar 
Blasche, S., Kim, Y., Oliveira, A. P. & Patil, K. R. Model microbial communities for ecosystems biology. Curr. Opin. Syst. Biol. 6, 51–57 (2017).Article 

Google Scholar 
Levy, S. F. et al. Quantitative evolutionary dynamics using high-resolution lineage tracking. Nature 519, 181–186 (2015).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Venkataram, S. et al. Development of a comprehensive genotype-to-fitness map of adaptation-driving mutations in yeast. Cell 166, 1585–1596.e22 (2016).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Jones, E. I., Bronstein, J. L. & Ferrière, R. The fundamental role of competition in the ecology and evolution of mutualisms. Ann. N. Y. Acad. Sci. 1256, 66–88 (2012).Article 
PubMed 

Google Scholar 
Boyer, S., Hérissant, L. & Sherlock, G. Adaptation is influenced by the complexity of environmental change during evolution in a dynamic environment. PLoS Genet. 17, e1009314 (2021).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Blundell, J. R. et al. The dynamics of adaptive genetic diversity during the early stages of clonal evolution. Nat. Ecol. Evol. 3, 293–301 (2019).Article 
PubMed 

Google Scholar 
Good, B. H., Rouzine, I. M., Balick, D. J., Hallatschek, O. & Desai, M. M. Distribution of fixed beneficial mutations and the rate of adaptation in asexual populations. Proc. Natl Acad. Sci. USA 109, 4950–4955 (2012).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Good, B. H., Martis, S. & Hallatschek, O. Adaptation limits ecological diversification and promotes ecological tinkering during the competition for substitutable resources. Proc. Natl Acad. Sci. USA 115, E10407–E10416 (2018).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Zhu, Y. O., Siegal, M. L., Hall, D. W. & Petrov, D. A. Precise estimates of mutation rate and spectrum in yeast. Proc. Natl Acad. Sci. USA 111, E2310–E2318 (2014).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Dunham, M. J. et al. Characteristic genome rearrangements in experimental evolution of Saccharomyces cerevisiae. Proc. Natl Acad. Sci. USA 99, 16144–16149 (2002).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Yona, A. H. et al. Chromosomal duplication is a transient evolutionary solution to stress. Proc. Natl Acad. Sci. USA 109, 21010–21015 (2012).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Sunshine, A. B. et al. The fitness consequences of aneuploidy are driven by condition-dependent gene effects. PLoS Biol. 13, e1002155 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Gerrish, P. J. & Lenski, R. E. The fate of competing beneficial mutations in an asexual population. Genetica 102-103, 127–144 (1998).Article 
CAS 
PubMed 

Google Scholar 
Desai, M. M. & Fisher, D. S. Beneficial mutation–selection balance and the effect of linkage on positive selection. Genetics 176, 1759–1798 (2007).Article 
PubMed 
PubMed Central 

Google Scholar 
Schiffels, S., Szöllosi, G. J., Mustonen, V. & Lässig, M. Emergent neutrality in adaptive asexual evolution. Genetics 189, 1361–1375 (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
Nguyen, Ba,A. N. et al. High-resolution lineage tracking reveals travelling wave of adaptation in laboratory yeast. Nature 575, 494–499 (2019).Article 

Google Scholar 
Foster, K. R., Shaulsky, G., Strassmann, J. E., Queller, D. C. & Thompson, C. R. L. Pleiotropy as a mechanism to stabilize cooperation. Nature 431, 693–696 (2004).Article 
CAS 
PubMed 

Google Scholar 
Sachs, J. L., Mueller, U. G., Wilcox, T. P. & Bull, J. J. The evolution of cooperation. Q. Rev. Biol. 79, 135–160 (2004).Article 
PubMed 

Google Scholar 
Harcombe, W. Novel cooperation experimentally evolved between species. Evolution 64, 2166–2172 (2010).PubMed 

Google Scholar 
Vasi, F., Travisano, M. & Lenski, R. E. Long-term experimental evolution in Escherichia coli. II. Changes in life-history traits during adaptation to a seasonal environment. Am. Nat. 144, 432–456 (1994).Article 

Google Scholar 
MacArthur, R. H. & Wilson, E. O. The Theory of Island Biogeography (Princeton Univ. Press, 2001).Reznick, D., Bryant, M. J. & Bashey, F. r- and K-selection revisited: the role of population regulation in life-history evolution. Ecology 83, 1509–1520 (2002).Article 

Google Scholar 
Mueller, L. D. & Ayala, F. J. Trade-off between r-selection and K-selection in Drosophila populations. Proc. Natl Acad. Sci. USA 78, 1303–1305 (1981).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Novak, M., Pfeiffer, T., Lenski, R. E., Sauer, U. & Bonhoeffer, S. Experimental tests for an evolutionary trade-off between growth rate and yield in E. coli. Am. Nat. 168, 242–251 (2006).Article 
PubMed 

Google Scholar 
Bachmann, H. et al. Availability of public goods shapes the evolution of competing metabolic strategies. Proc. Natl Acad. Sci. USA 110, 14302–14307 (2013).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Lipson, D. A. The complex relationship between microbial growth rate and yield and its implications for ecosystem processes. Front. Microbiol. 6, 615 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Orivel, J. et al. Trade-offs in an ant–plant–fungus mutualism. Proc. Biol. Sci. 284, 20161679 (2017).PubMed 
PubMed Central 

Google Scholar 
Fritts, R. K. et al. Enhanced nutrient uptake is sufficient to drive emergent cross-feeding between bacteria in a synthetic community. ISME J. 14, 2816–2828 (2020).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Wortel, M. T., Noor, E., Ferris, M., Bruggeman, F. J. & Liebermeister, W. Metabolic enzyme cost explains variable trade-offs between microbial growth rate and yield. PLoS Comput. Biol. 14, e1006010 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Cheng, C. et al. Laboratory evolution reveals a two-dimensional rate-yield tradeoff in microbial metabolism. PLoS Comput. Biol. 15, e1007066 (2019).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Luckinbill, L. S. r and K selection in experimental populations of Escherichia coli. Science 202, 1201–1203 (1978).Article 
CAS 
PubMed 

Google Scholar 
Oxman, E., Alon, U. & Dekel, E. Defined order of evolutionary adaptations: experimental evidence. Evolution 62, 1547–1554 (2008).Article 
PubMed 

Google Scholar 
Jasmin, J.-N., Dillon, M. M. & Zeyl, C. The yield of experimental yeast populations declines during selection. Proc. Biol. Sci. 279, 4382–4388 (2012).PubMed 
PubMed Central 

Google Scholar 
Laan, L., Koschwanez, J. H. & Murray, A. W. Evolutionary adaptation after crippling cell polarization follows reproducible trajectories. eLife 4, e09638 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Blount, Z. D., Lenski, R. E. & Losos, J. B. Contingency and determinism in evolution: replaying life’s tape. Science 362, eaam5979 (2018).Article 
PubMed 

Google Scholar 
Fukami, T. Historical contingency in community assembly: integrating niches, species pools, and priority effects. Annu. Rev. Ecol. Evol. Syst. 46, 1–23 (2015).Article 

Google Scholar 
Rainey, P. B. & Travisano, M. Adaptive radiation in a heterogeneous environment. Nature 394, 69–72 (1998).Article 
CAS 
PubMed 

Google Scholar 
Meyer, J. R. et al. Repeatability and contingency in the evolution of a key innovation in phage lambda. Science 335, 428–432 (2012).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Herron, M. D. & Doebeli, M. Parallel evolutionary dynamics of adaptive diversification in Escherichia coli. PLoS Biol. 11, e1001490 (2013).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Hillesland, K. L. et al. Erosion of functional independence early in the evolution of a microbial mutualism. Proc. Natl Acad. Sci. USA 111, 14822–14827 (2014).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Meroz, N., Tovi, N., Sorokin, Y. & Friedman, J. Community composition of microbial microcosms follows simple assembly rules at evolutionary timescales. Nat. Commun. 12, 2891 (2021).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
MacLean, R. C. The tragedy of the commons in microbial populations: insights from theoretical, comparative and experimental studies. Heredity 100, 471–477 (2008).Article 
CAS 
PubMed 

Google Scholar 
Dunn, B. et al. Recurrent rearrangement during adaptive evolution in an interspecific yeast hybrid suggests a model for rapid introgression. PLoS Genet. 9, e1003366 (2013).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Barillot, E., Lacroix, B. & Cohen, D. Theoretical analysis of library screening using a N-dimensional pooling strategy. Nucleic Acids Res. 19, 6241–6247 (1991).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Baym, M., Shaket, L., Anzai, I. A., Adesina, O. & Barstow, B. Rapid construction of a whole-genome transposon insertion collection for Shewanella oneidensis by Knockout Sudoku. Nat. Commun. 7, 13270 (2016).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Baym, M. et al. Inexpensive multiplexed library preparation for megabase-sized genomes. PLoS ONE 10, e0128036 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Venkataram, S., Kuo, H., Hom, E., Kryazhimskiy, S. Early adaptation in a microbial community is dominated by mutualism-enhancing mutations. Dryad https://doi.org/10.6076/D14K5X (2022). More

Plant species and pollinatorsMedicago sativa L. (Fabaceae), also called alfalfa or lucerne, is a perennial legume with flowers arranged in a cluster or raceme. It is a self-compatible plant with fairly high outcrossing rate (5.3–30%)46, and it requires insect visits for seed production47. No plant material was collected for this study. Honey bees, Apis mellifera, and alfalfa leafcutting bees, Megachile rotundata, are used as managed pollinators in alfalfa seed-production fields in the USA while bumble bees are commonly used in alfalfa breeding47.Experimental design and pollinator observationsFive 11 m × 11 m patches of M. sativa plants were set up in an east–west linear arrangement at the West Madison Agricultural Research Station in Madison, Wisconsin, USA. Within each patch, we transplanted 169 young plants grown from seeds in the greenhouse, each placed 90 cm apart. These plants grew and, at flowering, a plant had an average of 30.65 ± 16.4 stems per plant, with 4.93 ± 3.41 racemes per stem, and 7.53 ± 2.44 open flowers per raceme.A honey bee hive was placed approximately 100 m from the patches and a bumble bee hive was set up at the center of the southern edge of the patches. For leafcutting bees, a 60 × 30 × 7.6 cm bee board was set up in each of two boxes placed 1/3 and 2/3 along the southern edge of the patches and a half gallon of bees was released at periodic intervals throughout the alfalfa flowering season.Over two consecutive summers, observers followed bees foraging in the alfalfa patches, marked each raceme visited in succession within a foraging bout with a numbered clip, and recorded the number of flowers visited per raceme. After a bee had left a patch, observers went back to the marked racemes and measured the distance and direction traveled between consecutive racemes. Directions were recorded as one of the cardinal directions: North (N), South (S), East (E) or West (W), or inter-cardinal directions: Northeast (NE), Southeast (SE), Northwest (NW) and Southwest (SW). The frequency distributions of distances and directions traveled between two successive racemes are presented for each bee species each year in Figs. 1 (distances) and 2 (directions). The low pollinator abundance permitted observers to follow individual bees foraging in a patch. Little interference among bee species was observed in the patches.Figure 1Frequency distributions for distances traveled between consecutive racemes (cm) for each bee species each year.Full size imageFigure 2Frequency distributions of directions traveled between consecutive racemes for each bee species each year.Full size imageModel for the distance traveled between consecutive racemesWe first determined whether a statistical model best described the distance traveled between consecutive racemes (Modeled Distance), and examined whether the model differed among bee species. We used mixed effect linear models (proc Mixed in SAS 9.3)48 to identify the model that best described the distance traveled by pollinators between consecutive racemes. The model included loge distance as a linear function of loge flower number and bee species as fixed effects. The distance traveled between consecutive racemes and the number of flowers visited per raceme were log transformed prior to analyses in order to improve the models’ residuals. In addition, we included patch and foraging bout as random effects in the model. A foraging bout includes the racemes visited in succession from the time a bee is spotted in a patch to the time it leaves that patch. We used foraging bout instead of individual bee as the random effect because bees were not individually marked in this study. Moreover, to take into consideration the potential correlation between successive observations within a foraging bout, we added clip to the model. Clip 1 represents the first and second racemes visited in the foraging bout; clip 2, the second and third, and so on. Clip was added to the model either as a random effect or as a repeated measure with an AR(1) structure. The combination of random clip and random foraging bout creates a model that is sometimes called the “compound symmetry” model. The AR(1) structure represents correlations that decline exponentially as the gap between measurements increases such that measurements closer together in time are more strongly correlated than measurements further apart. Because we expected bees to visit flowers at close proximity when resources are abundant, we chose this correlation structure as a good potential descriptor of the way distances might be correlated within foraging bouts. We started with a full model which included loge flower number, bee species, patch, foraging bout, and clip either as a random effect or as a repeated measure with an AR(1) structure. We then removed variables and compared models by inspecting AIC values and the p values for each term in the model. We considered both low AIC and statistically significant (p  More

LightGBMBefore explaining LightGBM23, it is necessary to introduce XGBoost24, which is also based on the gradient boosting decision tree (GBDT) algorithm30. XGBoost integrates multiple classification and regression trees (CART) to compensate for the lack of prediction accuracy of a single CART. It is an improved boosting algorithm based on GBDT, which is popular due to its high processing speed, high regression accuracy and ability to process large-scale data31. However, XGBoost uses a presorted algorithm to find data segmentation points, which takes up considerable memory in the calculation and seriously affects cache optimization.LightGBM is improved based on XGBoost. It uses a histogram algorithm to find the best data segmentation point, which occupies less memory and has a lower complexity of data segmentation. The flow of the histogram algorithm to find the optimal segmentation point is shown in Fig. 3.Figure 3Histogram algorithm.Full size imageMoreover, LightGBM abandons the levelwise decision tree growth strategy used by most GBDT tools and uses the leafwise algorithm with depth limitations. This leaf-by-leaf growth strategy can reduce more errors and obtain better accuracy. Decision trees in boosting algorithms may grow too deep while training, leading to model overfitting. Therefore, LightGBM adds a maximum depth limit to the leafwise growth strategy to prevent this from happening and maintains its high computational efficiency. To summarize, LightGBM can be better and faster used in industrial practice and is also very suitable as the base model in our tide level prediction task. The layer-by-layer growth strategy and leaf-by-leaf growth strategy are shown in Fig. 4.Figure 4Two GBDT growth strategies.Full size imageCNN-BiGRUConvolutional neural networkA convolutional neural network (CNN) is a deep feedforward neural network with the characteristics of local connection and weight sharing. It was first used in the field of computer vision and achieved great success32,33. In recent years, CNNs have also been widely used in time series processing. For example, Bai et al.34 proposed a temporal convolutional network (TCN) based on a convolutional neural network and residual connections, which is not worse than recurrent neural networks such as LSTM in some time series analysis tasks. At present, a convolutional neural network is generally composed of convolution layers, pooling layers and a fully connected layer. Its network structure is shown in Fig. 5. The pooling layer is usually added after the convolution layers. The maximum pooling layer can retain the strong features in the data after the convolution operation, eliminate the weak features to reduce the number of parameters in a network and avoid overfitting of the model.Figure 5Schematic diagram of a convolutional neural network.Full size imageBidirectional GRUIn previous attempts at tide level prediction by scholars, bidirectional long short-term memory networks35 have achieved good prediction results. However, in our subsequent experiments, the bidirectional gated recurrent unit achieved higher prediction accuracy than BiLSTM, so we used the BiGRU network for subsequent prediction tasks.The GRU network36 adds a gating mechanism to control information updating in a recurrent neural network. Different from the mechanism in LSTM, GRU consists of only two gates called the update gate ({z}_{t}) and the reset door ({r}_{t}).The recurrent unit structure of the GRU network is shown in Fig. 6.Figure 6Recurrent unit structure of the GRU network.Full size imageEach unit of GRU is calculated as follows:$${z}_{t}= sigma ({W}_{z}{x}_{t}+{U}_{z}{h}_{t-1}+{b}_{z})$$
(7)
$${r}_{t}= sigma ({W}_{r}{x}_{t}+{U}_{r}{h}_{t-1}+{b}_{r})$$
(8)
$${widetilde{h}}_{t}=tanh({W}_{h}{x}_{t}+{U}_{h}left({r}_{t}odot {h}_{t-1}right)+{b}_{h})$$
(9)
$${h}_{t}={z}_{t}odot {h}_{t-1}+left(1-{z}_{t}right)odot {widetilde{h}}_{t}$$
(10)
In the above formula, ({z}_{t}) represents the update gate, which controls how much information is retained from the previous state ({h}_{t-1}) (without nonlinear transformation) when calculating the current state ({h}_{t}). Meanwhile, it also controls how much information will be accepted by ({h}_{t}) from the candidate states ({widetilde{h}}_{t}). ({r}_{t}) represents the reset gate, which is used to ensure whether the calculation of the candidate state ({widetilde{h}}_{t}) depends on the previous state ({h}_{t-1}). (upsigma ) is the standard sigmoid activation function; (tanh(cdot )) is the hyperbolic tangent activation function; and (odot ) indicates the Hadamard product. The weight matrices of the update gate, reset gate, and ({widetilde{h}}_{t}) calculation layer are expressed as ({W}_{z},{W}_{r},{W}_{h}); the coefficient matrices of the update gate, reset gate, and ({widetilde{h}}_{t}) calculation layer are expressed as ({U}_{z},{U}_{r},{U}_{h}); and the offset vectors of the update gate, reset gate, and ({widetilde{h}}_{t}) calculation layer are expressed as ({b}_{z},{b}_{r},{b}_{h}).A bidirectional gated recurrent unit network37 is a combination of two GRUs whose information propagating directions are reversed, and it has independent parameters in each, which makes it able to fit both forward and backward data at first and then join up the results from two directions. BiGRU can capture sequence patterns that may be ignored by unidirectional GRU. The structure of BiGRU is shown in Fig. 7.Figure 7The structure of BiGRU.Full size imageTaking the BiGRU’s forward hidden state vector at time (t) as ({h}_{t}^{(1)}) and taking the BiGRU’s backward hidden state vector at time (t) as ({h}_{t}^{(2)}), (upsigma ) indicates the standard sigmoid activation function, and (oplus ) indicates a vector splicing operation. We can calculate the output ({y}_{t}) of a BiGRU network as follows:$${h}_{t}={h}_{t}^{(1)}oplus {h}_{t}^{(2)}$$
(11)
$${y}_{t}=sigma ({h}_{t} )$$
(12)
CNN-BiGRU prediction modelBecause CNN has significant advantages in extracting useful features from a picture or a sequence and BiGRU is good at processing time series, we combine CNN and BiGRU to build the CNN-BiGRU model. The model can be mainly divided into an input layer, a convolution layer, a BiGRU network layer, a dropout layer, a fully connected layer and an output layer. The CNN layer and BiGRU layer are the core structures of the model. The function of the dropout layer is to avoid model overfitting. The CNN layer consists of two one-dimensional convolution (Conv1D) layers and a one-dimensional maximum pooling (MaxPooling1D) layer. The input of BiGRU is the output sequence of the CNN layer, and the BiGRU network is set as a one-hidden-layer structure. The structure of the CNN-BiGRU combination model is shown in Fig. 8.Figure 8The structure of CNN-BiGRU.Full size imageVariable weight combination modelWhen we analyze and predict relatively stationary tide level time series, LightGBM can perform well. However, due to environmental factors such as air pressure, wind force and terrain in reality, most tide level observation sequences are sometimes not relatively stationary, which requires that our tide level prediction model be reasonably able to “extrapolate” based on the sample observations, that is, be capable of generating values that are not in the sample. LightGBM is a tree-based model, which leads to our prediction results being between the maximum and minimum values of sequences. Therefore, LightGBM will not be able to accurately predict the situation or tidal change trend that did not appear in previous observations. However, the CNN-BiGRU model, which is a kind of neural network, has no such problem in theory and will be able to find the trend information that may be hidden in the tide level series. Therefore, we consider providing an appropriate weight for a single base model to build a combination model to improve the accuracy of the tide level prediction task.Principle of the residual weight combination model and improved variable weight combination modelTo improve the prediction accuracy of the combination model, a simple and effective idea is to determine the base models’ weights in the combination model according to the error between the prediction value and the real value. This method is also called the residual weight method, and its calculation formulas for determining the weights are:$$gleft({x}_{t}right)= sum_{i=1}^{m}{omega }_{i}left(t-1right){f}_{i}({x}_{t})$$
(13)
$${omega }_{i}left(t-1right)=frac{frac{1}{overline{{varphi }_{i}}left(t-1right)}}{sum_{i=1}^{m}frac{1}{overline{{varphi }_{i}}left(t-1right)}}$$
(14)
$$sum_{i=1}^{m}{omega }_{i}left(t-1right)=1,{omega }_{i}left(t-1right)ge 0$$
(15)

where ({omega }_{i}left(t-1right)) denotes the weight of the (i) th model at the moment (t-1), ({f}_{i}left({x}_{t}right)) denotes the prediction value of the (i) th model at the moment (t), (gleft({x}_{t}right)) denotes the prediction value of the combination model at the moment (t), and (overline{{varphi }_{i}}left(t-1right)) is the square sum of the predictive errors of the (i) th model at the moment (t-1).Our LightGBM-CNN-BiGRU (combination model) is based on the improved residual weight method. We call it the variable weight combination model. We use the weights calculated by formula (9) and formula (11) to calculate a series of new weights. The new weights from formula (11) will take the residual weight changes in (d) time steps into consideration by averaging the old weights in (d) time steps to improve the stability of the residual weight method.$${omega }_{j}left(tright)=frac{1}{d}sum_{k=1}^{d}{omega }_{i}left(t-kright)left(d=4right)$$
(16)
After obtaining a series of weights through formula (9) and formula (11), we take the absolute value of the error between the prediction value and the true value of each combination model at the moment of (t) as ({delta }_{i,t}) and ({delta }_{j,t}), respectively:$${delta }_{i,t}=mid sum_{i=1}^{m}{omega }_{i}left(tright){f}_{i}left({x}_{t}right)-{y}_{t}mid $$
(17)
$${delta }_{j,t}=mid sum_{i=1}^{m}{omega }_{j}left(tright){f}_{i}left({x}_{t}right)-{y}_{t}mid $$
(18)
Comparing ({delta }_{i,t}) and ({delta }_{j,t}), if ({delta }_{i,t} >{delta }_{j,t}), the combination model uses the new weight ({omega }_{j}left(tright)) in place of the original weight ({omega }_{i}left(tright)). Otherwise, the weight of the combination model remains unchanged.Parameter optimization of the combination modelBecause the LightGBM-CNN-BiGRU (combination model) is a variable weight combination of the prediction results from two base models, the performance of the combination model can be directly improved by separately optimizing the super parameters of the two base models. We mainly use the grid search algorithm and K-fold cross validation method to optimize the parameters. The grid search algorithm is a method to improve the performance of a certain model by iterating over a given set of parameters. With the help of the K-fold cross validation method, we can calculate the performance score of the LightGBM model on the training set and easily optimize its superparameters. The final parameters of the LightGBM model are set to num_leaves = 26, learning_rate = 0.05, and n_estimators = 46.For the CNN-BiGRU network, we mainly improve the prediction accuracy of the model by adjusting the size and number of hidden layers in the BiGRU structure and prevent the model from overfitting by changing the dropout ratio and tracking the validation loss of the network while training.The LightGBM and CNN-BiGRU variable weight combination modelThe workflow of our tide level prediction model is shown in Fig. 9. It mainly includes data preprocessing; training, optimization and prediction of the base models; construction of a variable weight combination prediction model; and evaluation and analysis of the combination model’s performance.

(1)

Data preprocessing: The quality of the data directly determines the upper limit of the prediction and generalization ability of a certain machine learning model. Standard, clean and continuous data are conducive to model training. The data used in this study are from the Irish National Tide Gauge Network, and all of them are subject to quality control. We filled in a small number of missing values and normalized the data to speed up the model training.

(2)

Construction and optimization of base models: We divide the dataset into a training set, a validation set and a test set according to the proportion of 7:1:2 and train the LightGBM model and CNN-BiGRU model with data on the training set. We optimize the parameters and monitor whether the model has been overfitted by tracking the validation loss of the network while training. Finally, we put the data into two base models for training and then obtain the prediction results of a single base model.

(3)

Construction of the variable weight combination model. Based on the prediction results of two single base models obtained in step (2), we calculate the weight of each base model according to the principle of the improved variable weight combination method and then obtain the prediction results of the variable weight combination model.

(4)

Model evaluation and analysis: According to the indexes of the model evaluation, the variable weight combination model is compared with other basic models to analyze its prediction performance after being improved.

Figure 9Prediction flow of the LightGBM-CNN-BiGRU variable weight combination model.Full size image More

Snow AJD, Burchill L, Sharma M, Davies GJ, Williams SJ. Sulfoglycolysis: Catabolic pathways for metabolism of sulfoquinovose. Chem Soc Rev. 2021;50:13628–45.Article 
CAS 

Google Scholar 
Van Mooy BAS, Rocap G, Fredricks HF, Evans CT, Devol AH. Sulfolipids dramatically decrease phosphorus demand by picocyanobacteria in oligotrophic marine environments. Proc Natl Acad Sci USA 2006;103:8607–12.Article 

Google Scholar 
Wu J, Sunda W, Boyle EA, Karl DM. Phosphate depletion in the western North Atlantic. Ocean Sci 2000;289:759–62.CAS 

Google Scholar 
Goddard-Borger ED, Williams SJ. Sulfoquinovose in the biosphere: occurrence, metabolism and functions. Biochem J. 2017;474:827–49.Article 
CAS 

Google Scholar 
Harwood JL, Nicholls RG. The plant sulpholipid- a major component of the sulphur cycle. Biochem Soc Trans. 1979;7:440–7.Article 
CAS 

Google Scholar 
Moran MA, Durham BP. Sulfur metabolites in the pelagic ocean. Nat Rev Microbiol. 2019;17:665–78.Article 
CAS 

Google Scholar 
Tang K. Chemical diversity and biochemical transformation of biogenic organic sulfur in the ocean. Front Mar Sci. 2020;7:68.Article 

Google Scholar 
Denger K, Weiss M, Felux AK, Schneider A, Mayer C, Spiteller D, et al. Sulphoglycolysis in Escherichia coli K-12 closes a gap in the biogeochemical sulphur cycle. Nature 2014;507:114–7.Article 
CAS 

Google Scholar 
Hanson BT, Kits KD, Loffler J, Burrichter AG, Fiedler A, Denger K, et al. Sulfoquinovose is a select nutrient of prominent bacteria and a source of hydrogen sulfide in the human gut. ISME J. 2021;15:2779–91.Article 
CAS 

Google Scholar 
Strickland TC, Fitzgerald JW. Mineralization of sulfur in sulfoquinovose by forest soils. Soil Biol Biochem. 1983;15:347–9.Article 
CAS 

Google Scholar 
Felux AK, Spiteller D, Klebensberger J, Schleheck D. Entner-Doudoroff pathway for sulfoquinovose degradation in Pseudomonas putida SQ1. Proc Natl Acad Sci USA 2015;112:E4298–E305.Article 
CAS 

Google Scholar 
Frommeyer B, Fiedler AW, Oehler SR, Hanson BT, Loy A, Franchini P, et al. Environmental and intestinal phylum Firmicutes bacteria metabolize the plant sugar sulfoquinovose via a 6-deoxy-6-sulfofructose transaldolase pathway. Iscience. 2020;23:101510.Article 
CAS 

Google Scholar 
Roy AB, Hewlins MJE, Ellis AJ, Harwood JL, White GF. Glycolytic breakdown of sulfoquinovose in bacteria: A missing link in the sulfur cycle. Appl Environ Microbiol. 2003;69:6434–41.Article 
CAS 

Google Scholar 
Liu J, Wei Y, Ma K, An J, Liu X, Liu Y, et al. Mechanistically diverse pathways for sulfoquinovose degradation in bacteria. ACS Catal. 2021;11:14740–50.Article 
CAS 

Google Scholar 
Zhang S, Li Z, Yan Y, Zhang C, Li J, Zhao B. Bacillus urumqiensis sp. nov., a moderately haloalkaliphilic bacterium isolated from a salt lake. Int J Syst Evol Microbiol. 2016;66:2305–12.Article 
CAS 

Google Scholar 
Durham BP, Sharma S, Luo H, Smith CB, Amin SA, Bender SJ, et al. Cryptic carbon and sulfur cycling between surface ocean plankton. Proc Natl Acad Sci USA 2015;112:453–7.Article 
CAS 

Google Scholar 
Chen X, Liu L, Gao X, Dai X, Han Y, Chen Q, et al. Metabolism of chiral sulfonate compound 2,3-dihydroxypropane-1-sulfo-nate (DHPS) by Roseobacter bacteria in marine environment. Environ Int. 2021;157:106829.Article 
CAS 

Google Scholar 
Liu J, Wei Y, Lin L, Teng L, Yin J, Lu Q, et al. Two Radical-dependent mechanisms for anaerobic degradation of the globally abundant Organosulfur Compound Dihydroxypropanesulfonate. Proc Natl Acad Sci USA 2020;117:15599.Article 
CAS 

Google Scholar 
Xing M, Wei Y, Zhou Y, Zhang J, Lin L, Hu Y, et al. Radical-mediated C-S bond cleavage in C2 sulfonate degradation by anaerobic bacteria. Nat Commun. 2019;10:1609.Article 

Google Scholar 
Sharma M, Lingford JP, Petricevic M, Snow AJD, Zhang Y, Jarva MA, et al. Oxidative desulfurization pathway for complete catabolism of sulfoquinovose by bacteria. Proc Natl Acad Sci USA 2022;119:e2116022119.Article 
CAS 

Google Scholar 
Scholz SS, Serif M, Schleheck D, Sayer MDJ, Cook AM, Kupper FC. Sulfoquinovose metabolism in marine algae. Bot Mar. 2021;64:301–12.Article 
CAS 

Google Scholar 
Abayakoon P, Epa R, Petricevic M, Bengt C, Mui JWY, van der Peet PL, et al. Comprehensive synthesis of substrates, intermediates, and products of the sulfoglycolytic Embden-Meyerhoff-Parnas pathway. J Org Chem. 2019;84:2901–10.Article 
CAS 

Google Scholar 
Denger K, Smits THM, Cook AM. L-Cysteate sulpho-lyase, a widespread pyridoxal 5 ‘-phosphate-coupled desulphonative enzyme purified from Silicibacter pomeroyi DSS-3. Biochem J. 2006;394:657–64.Article 
CAS 

Google Scholar 
Guillard RRL. Culture of Phytoplankton for Feeding Marine Invertebrates. Smith WL, Chanley MH, (eds): Springer US; 1975. Boston, MA. pp 29–60.Moore LR, Coe A, Zinser ER, Saito MA, Sullivan MB, Lindell D, et al. Culturing the marine cyanobacterium Prochlorococcus. Limnol Oceanogr Methods. 2007;5:353–62.Article 
CAS 

Google Scholar 
Waterbury J, Watson S, Valois F, Franks D. Biological and ecological characterization of the marine unicellular cyanobacterium Synechococcus. Platt T, Li WKW, (eds). Department of Fisheries and Oceans, Ottawa 1986. pp 71–120.Olenina I, Hajdu S, Edler L, Andersson A, Wasmund N, Busch S, et al. Biovolumes and size-classes of phytoplankton in the Baltic Sea. HELCOM Balt Sea Environ Proc. 2006;106:144.
Google Scholar 
Zheng Q, Wang Y, Lu J, Lin W, Chen F, Jiao N. Metagenomic and metaproteomic insights into photoautotrophic and heterotrophic interactions in a Synechococcus culture. mbio 2020;11:e03261–19.Article 
CAS 

Google Scholar 
Partensky F, Hess WR, Vaulot D. Prochlorococcus, a marine photosynthetic prokaryote of global significance. Microbiol Mol Biol Rev. 1999;63:106–27.Article 
CAS 

Google Scholar 
Han Y, Zhang M, Chen X, Zhai W, Tan E, Tang K. Transcriptomic evidences for microbial carbon and nitrogen cycles in the deoxygenated seawaters of Bohai Sea. Environ Int. 2022;158:106889.Article 
CAS 

Google Scholar 
Li WKW. Primary production of prochlorophytes, cyanobacteria, and eukaryotic ultraphytoplankton – measurements from flow cytometric sorting. Limnol Oceanogr. 1994;39:169–75.Article 
CAS 

Google Scholar 
Denger K, Ruff A, Rein U, Cook AM. Sulphoacetaldehyde sulpho-lyase (EC 4.4.1.12) from Desulfonispora thiosulfatigenes: purification, properties and primary sequence. Biochem J. 2001;357:581–6.Article 
CAS 

Google Scholar 
Ismail R, Lee HY, Mahyudin NA, Abu, Bakar F. Linearity study on detection and quantification limits for the determination of avermectins using linear regression. J Food Drug Anal. 2014;22:407–12.Article 
CAS 

Google Scholar 
Klemetsen T, Raknes IA, Fu J, Agafonov A, Balasundaram SV, Tartari G, et al. The MAR databases: development and implementation of databases specific for marine metagenomics. Nucleic Acids Res. 2018;46:D692–D9.Article 
CAS 

Google Scholar 
Suzek BE, Huang H, McGarvey P, Mazumder R, Wu CH. UniRef: comprehensive and non-redundant UniProt reference clusters. Bioinformatics 2007;23:1282–8.Article 
CAS 

Google Scholar 
Rozewicki J, Li S, Amada KM, Standley DM, Katoh K. MAFFT-DASH: Integrated protein sequence and structural alignment. Nucleic Acids Res. 2019;47:W5–W10.CAS 

Google Scholar 
Schuller DJ, Reisch CR, Moran MA, Whitman WB, Lanzilotta WN. Structures of dimethylsulfoniopropionate-dependent demethylase from the marine organism Pelagabacter ubique. Protein Sci. 2012;21:289–98.Article 
CAS 

Google Scholar 
Bharath SR, Bisht S, Harijan RK, Savithri HS, Murthy MR. Structural and mutational studies on substrate specificity and catalysis of Salmonella typhimurium D-cysteine desulfhydrase. PLoS One. 2012;7:e36267.Article 
CAS 

Google Scholar 
Chartron J, Carroll KS, Shiau C, Gao H, Leary JA, Bertozzi CR, et al. Substrate Recognition, Protein Dynamics, and Iron-Sulfur Cluster in Pseudomonas aeruginosa Adenosine 5′-Phosphosulfate Reductase. J Mol Biol. 2006;364:152–69.Article 
CAS 

Google Scholar 
Davis KM, Altmyer M, Martinie RJ, Schaperdoth I, Krebs C, Bollinger JM Jr, et al. Structure of a Ferryl Mimic in the Archetypal Iron(II)- and 2-(Oxo)-glutarate-Dependent Dioxygenase, TauD. Biochemistry 2019;58:4218–23.Article 
CAS 

Google Scholar 
Mirdita M, Schütze K, Moriwaki Y, Heo L, Ovchinnikov S, Steinegger M. ColabFold: Making protein folding accessible to all. Nat Methods. 2022;19:679–82.Article 
CAS 

Google Scholar 
Jumper J, Evans R, Pritzel A, Green T, Figurnov M, Ronneberger O, et al. Highly accurate protein structure prediction with AlphaFold. Nature 2021;596:583–9.Article 
CAS 

Google Scholar 
Zhang C, Shine M, Pyle AM, Zhang Y. US-align: universal structure alignments of proteins, nucleic acids, and macromolecular complexes. Nat Methods. 2022;19:1109–15.Article 
CAS 

Google Scholar 
Xu J, Zhang Y. How significant is a protein structure similarity with TM-score = 0.5? Bioinformatics 2010;26:889–95.Article 
CAS 

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

Google Scholar 
Villar E, Vannier T, Vernette C, Lescot M, Cuenca M, Alexandre A, et al. The Ocean Gene Atlas: exploring the biogeography of plankton genes online. Nucleic Acids Res. 2018;46:W289–W95.Article 
CAS 

Google Scholar 
Vernette C, Henry N, Lecubin J, de Vargas C, Hingamp P, Lescot M. The Ocean barcode atlas: A web service to explore the biodiversity and biogeography of marine organisms. Mol Ecol Resour. 2021;21:1347–58.Article 
CAS 

Google Scholar 
Paoli L, Ruscheweyh H-J, Forneris CC, Hubrich F, Kautsar S, Bhushan A, et al. Biosynthetic potential of the global ocean microbiome. Nature 2022;607:111–8.Article 
CAS 

Google Scholar 
Sunagawa S, Acinas SG, Bork P, Bowler C, Acinas SG, Babin M, et al. Tara Oceans: towards global ocean ecosystems biology. Nat Rev Microbiol. 2020;18:428–45.Article 
CAS 

Google Scholar 
Acinas SG, Sánchez P, Salazar G, Cornejo-Castillo FM, Sebastián M, Logares R, et al. Deep ocean metagenomes provide insight into the metabolic architecture of bathypelagic microbial communities. Commun Biol. 2021;4:604.Article 
CAS 

Google Scholar 
Biller SJ, Berube PM, Dooley K, Williams M, Satinsky BM, Hackl T, et al. Marine microbial metagenomes sampled across space and time. Sci Data. 2018;5:180176.Article 
CAS 

Google Scholar 
Pachiadaki MG, Brown JM, Brown J, Bezuidt O, Berube PM, Biller SJ, et al. Charting the Complexity of the Marine Microbiome through Single-Cell Genomics. Cell 2019;179:1623–35.Article 
CAS 

Google Scholar 
Delmont TO, Quince C, Shaiber A, Esen ÖC, Lee STM, Rappé MS, et al. Nitrogen-fixing populations of Planctomycetes and Proteobacteria are abundant in surface ocean metagenomes. Nat Microbiol. 2018;3:804–13.Article 
CAS 

Google Scholar 
Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: Molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30:2725–9.Article 
CAS 

Google Scholar 
Subramanian B, Gao S, Lercher MJ, Hu S, Chen W-H. Evolview v3: A webserver for visualization, annotation, and management of phylogenetic trees. Nucleic Acids Res. 2019;47:W270–W5.Article 
CAS 

Google Scholar 
Xing M, Wei Y, Zhou Y, Zhang J, Lin L, Hu Y, et al. Radical-mediated C-S bond cleavage in C2 sulfonate degradation by anaerobic bacteria. Nat Commun. 2019;10:1609.Article 

Google Scholar 
Biebl H, Allgaier M, Tindall BJ, Koblizek M, Lunsdorf H, Pukall R, et al. Dinoroseobacter shibae gen. nov., sp nov., a new aerobic phototrophic bacterium isolated from dinoflagellates. Int J Syst Evol Microbiol. 2005;55:1089–96.Article 
CAS 

Google Scholar 
Fu H, Uchimiya M, Gore J, Moran MA. Ecological drivers of bacterial community assembly in synthetic phycospheres. Proc Natl Acad Sci USA 2020;117:3656–62.Article 
CAS 

Google Scholar 
Chen I-MA, Chu K, Palaniappan K, Ratner A, Huang J, Huntemann M, et al. The IMG/M data management and analysis system v.7: content updates and new features. Nucleic Acids Res. 2022. https://doi.org/10.1093/nar/gkac976.Shiba T. Roseobacter litoralis gen. nov., sp. nov., and Roseobacter denitrificans sp. nov., aerobic pink-pigmented bacteria which contain bacteriochlorophyll a. Syst Appl Microbiol. 1991;14:140–5.Article 

Google Scholar 
Kopriva S, Calderwood A, Weckopp SC, Koprivova A. Plant sulfur and big data. Plant Sci. 2015;241:1–10.Article 
CAS 

Google Scholar 
Simon J, Kroneck PMH. Microbial sulfite respiration. Adv Micro Physiol. 2013;62:45–117.Article 
CAS 

Google Scholar 
Gonzalez JM, Covert JS, Whitman WB, Henriksen JR, Mayer F, Scharf B, et al. Silicibacter pomeroyi sp nov and Roseovarius nubinhibens sp nov., dimethylsulfoniopropionate-demethylating bacteria from marine environments. Int J Syst Evol Microbiol. 2003;53:1261–9.Article 
CAS 

Google Scholar 
Liang KYH, Orata FD, Boucher YF, Case RJ. Roseobacters in a sea of poly- and paraphyly: whole genome-based taxonomy of the family Rhodobacteraceae and the proposal for the split of the “Roseobacter clade” into a novel family, Roseobacteraceae fam. nov. Front Microbiol. 2021;12:683109.Article 

Google Scholar 
Howard EC, Sun S, Biers EJ, Moran MA. Abundant and diverse bacteria involved in DMSP degradation in marine surface waters. Environ Microbiol. 2008;10:2397–410.Article 
CAS 

Google Scholar 
Howard EC, Henriksen JR, Buchan A, Reisch CR, Buergmann H, Welsh R, et al. Bacterial taxa that limit sulfur flux from the ocean. Science. 2006;314:649–52.Article 
CAS 

Google Scholar 
Durham BP, Boysen AK, Carlson LT, Groussman RD, Heal KR, Cain KR, et al. Sulfonate-based networks between eukaryotic phytoplankton and heterotrophic bacteria in the surface ocean. Nat Microbiol. 2019;4:1706–15.Article 
CAS 

Google Scholar 
Smetacek V. Diatoms and the ocean carbon cycle. Protist 1999;150:25–32.Article 
CAS 

Google Scholar 
Stoecker DK, Lavrentyev PJ. Mixotrophic plankton in the polar seas: A pan-Arctic review. Front Mar Sci. 2018;5:292.Article 

Google Scholar 
Turner SM, Malin G, Liss PS, Harbour DS, Holligan PM. The seasonal-variation of dimethyl sulfide and dimethylsulfoniopropionate concentrations in nearshore waters. Limnol Oceanogr. 1988;33:364–75.Article 
CAS 

Google Scholar 
Belviso S, Kim S-K, Rassoulzadegan F, Krajka B, Nguyen BC, Mihalopoulos N, et al. Production of dimethylsulfonium propionate (DMSP) and dimethylsulfide (DMS) by a microbial food web. Limnol Oceanogr. 1990;35:1810–21.Article 
CAS 

Google Scholar 
Simo R, Pedros-Alio C, Malin G, Grimalt JO. Biological turnover of DMS, DMSP and DMSO in contrasting open-sea waters. Mar Ecol Prog Ser. 2000;203:1–11.Article 
CAS 

Google Scholar 
Flombaum P, Gallegos JL, Gordillo RA, Rincon J, Zabala LL, Jiao N, et al. Present and future global distributions of the marine Cyanobacteria Prochlorococcus and Synechococcus. Proc Natl Acad Sci USA 2013;110:9824–9.Article 
CAS 

Google Scholar 
Gasparovic B, Penezic A, Frka S, Kazazic S, Lampitt RS, Holguin FO, et al. Particulate sulfur-containing lipids: Production and cycling from the epipelagic to the abyssopelagic zone. Deep Sea Res Part I Oceanogr Res Pap. 2018;134:12–22.Article 
CAS 

Google Scholar 
Zhan P, Tang K, Chen X, Yu L. Complete genome sequence of Maribacter sp T28, a polysaccharide-degrading marine flavobacteria. J Biotechnol. 2017;259:1–5.Article 
CAS 

Google Scholar 
Van Mooy BAS, Fredricks HF. Bacterial and eukaryotic intact polar lipids in the eastern subtropical South Pacific: Water-column distribution, planktonic sources, and fatty acid composition. Geochim Cosmochim Acta. 2010;74:6499–516.Article 

Google Scholar 
Popendorf KJ, Tanaka T, Pujo-Pay M, Lagaria A, Courties C, Conan P, et al. Gradients in intact polar diacylglycerolipids across the Mediterranean Sea are related to phosphate availability. Biogeosciences 2011;8:3733–45.Article 
CAS 

Google Scholar  More

Restoration opportunity area and costsMangrove restoration programmes have a greater chance of being successful when implemented in areas where mangroves have previously grown15. These areas have either been subject to deforestation or degradation and may be under government management or private ownership. They are locations that have undergone forest conversion into other land uses, including aquaculture, crops or plantations and urban settlements. Land ownership status is an important factor to consider for determining the availability of land for mangrove restoration7. For example, a higher opportunity and priority would be given to unproductive aquaculture ponds located in the protected and production forest areas which are under government management or leasehold, rather than in areas with other land uses that may be under private ownership (Methods gives detailed forest land tenure classifications in Indonesia). Therefore, managing mangrove rehabilitation should consider factors that include land tenure status and land-cover type as well as biogeomorphology (for example, ensuring that the correct mangrove species are used in hydrologically suitable locations) across landscape scales.We calculated that ~193,367 ha of land may be feasible for implementation of mangrove rehabilitation programmes (Fig. 4). This conservative assessment suggests that the potential for restoration may be only 30% of the current mangrove rehabilitation area target (600,000 ha). Depending on the challenges and opportunities for each of the biogeomorphological categories of land use and the forest land status we considered (see Methods for detailed mapping methodology), we identified that 9% of the potential restorable area was categorized as being within the high opportunity scenario, 33% as medium and 58% as areas falling within the low opportunity scenario. Among these scenarios, ~75% of identified areas have non-protected forest status, implying a greater tenurial challenge to establishing a rehabilitation programme. We identified the five provinces that are among the top ranked of high potential for mangrove restoration in Indonesia, namely East Kalimantan (20% of national restoration potential area), North Kalimantan (20%), South Sumatra (12%), West Kalimantan (5%) and Riau provinces (5%) (Fig. 1c). All of these provinces, except South Sumatra, are among the areas already identified in the current mangrove rehabilitation programme by the BRGM as having high opportunity for rehabilitation4. At the subprovincial scale, we identified the top six regencies with restoration area opportunity >10,000 ha, namely Banyuasin, Bulungan, Tana Tidung, Paser, Berau and Nunukan (Supplementary Table 1). Mangroves across these regions were commonly deforested after 2010 and converted into aquaculture ponds despite being designated as protected forest areas (Supplementary Table 1).Fig. 4: The distribution of mangrove loss area (in hectares) between 2001 and 2020 in Indonesia.Also shown are mangrove loss proportions within different biogeomorphological typology, loss drivers (land-use types), forest land status and identified scenarios of restoration opportunity (low, medium and high).Full size imageConsidering that previous successful (85% survival rates) mangrove rehabilitation around the world has been achieved only at small landscape scales (10–400 ha) with costs varying between US$1,500 ha−1 and US$9,000 ha−1 (refs. 8,16), the large-scale mangrove rehabilitation ambition of Indonesia must be carefully planned. Rehabilitating ~200,000 ha of degraded mangroves will require between US$0.29 billion and US$1.74 billion. The 2021 annual government budget allocation for mangrove rehabilitation under BRGM alone is ~US$0.10 billion17, which is 66–94% lower than the estimated total required budget but with additional international investment18 there is potential for scalable mangrove rehabilitation success.Lessons learned from the past failuresIn Indonesia, unproductive aquaculture ponds have become targets for mangrove rehabilitation programmes (Supplementary Fig. 1). However, metrics of rehabilitation success in these settings reveal low survival rates of planted seedlings, highlighting an urgency to develop new strategies for mangrove rehabilitation and strategies to assess the effectiveness of ecosystem rehabilitation6. For example, a silviculture approach—nursery-based mangrove planting using Rhizophora species—has been adopted for mangrove restoration and management for a long time in Indonesia19. When seedlings are directly planted in unused ponds (Supplementary Fig. 1), dense monoculture plantations often form, which despite providing some ecosystem services (for example, carbon sequestration20) have limited biodiversity value21 and may be less resilient to stressors compared to a diverse assemblages of tree species22.Mangrove restoration projects have often suffered low success rates due to inadequate hydrological site assessments before revegetation23. For example, mangrove planting programmes initiated after the 2004 tsunami were focused on mono-species planting and on reporting the number of seedlings being planted in a given area24. These planting projects most often occurred on undisputed land, such as mudflats, which are inappropriate locations for long-term mangrove growth because of high inundation frequency, high water flow rates and hypersaline conditions that limit seedling establishment and survival24. Planting has also focused in mangrove areas where low canopy cover is observed. While some mangrove areas with low canopy cover may respond to plantings because they are degraded, many sites naturally support low canopy cover, reflecting suboptimal environmental conditions for growth of Rhizophora species, instead favouring growth of highly salt tolerant species such as Avicennia spp.24. Such failures in mangrove rehabilitation efforts, however, have been under-reported with more than 50% of rehabilitation studies not monitored over time (Supplementary Fig. 1).Alternative restoration approaches through repairing hydrology, including excavation and removal of pond walls and tidal gates, have also been introduced15, although this approach has been only practiced in Indonesia at limited scales, mostly in unused aquaculture ponds25. A comprehensive understanding of the opportunity for mangrove rehabilitation in Indonesia is largely unquantified. Additionally, with limited monitoring of mangrove rehabilitation projects, the effectiveness and functionality of mangrove rehabilitation in Indonesia remains largely unknown and therefore it remains challenging to assess rehabilitation effectiveness between approaches and locations in Indonesia. Yet such assessments provide important data to achieve the ambitious mangrove rehabilitation goals of Indonesia.Mangrove governance in IndonesiaMangrove conservation in Indonesia was formally adopted in 1990 (Extended Data Fig. 1 and Supplementary Table 2), when mangroves were designated as protected forests under Law 5/1990 and the Presidential Decree 32/1990. When the Asian tsunami hit Aceh province in 2004, the role of mangroves in wave attenuation and therefore minimizing disaster risks for coastal communities was recognized26. As a result, nearly 30,000 ha of damaged mangroves were rehabilitated to recover coastal resiliency through planting of nearly 24 million seedlings over 60 projects24. However, the success of these programmes was low due to a lack of planning, monitoring and critical supplemental actions24,27. Despite the failure of many mangrove rehabilitation projects post-tsunami, the implementation of the subsequent programmes have not fully adopted best-practice mangrove rehabilitation principles6,7,15,23. In 2007, similar approaches to mangrove rehabilitation and conservation were adopted at a larger, national scale under the Spatial Planning Law (Law 26/2007) and the Coastal Area and Small Islands Management Law (Law 27/2007).In 2012, the National Mangrove Management Strategy (STRANAS Mangrove) was first established and followed by the formalization of the National and Regional Mangrove Working Group whose task was to guide mangrove conservation and rehabilitation. Its main goal was to involve more stakeholders, including civil society organizations and subnational government bodies, in mangrove conservation and rehabilitation28. Until 2017, the technical regulation of strategy and performance indicators for mangrove management was implemented with targets set to rehabilitate 3.49 Mha of mangroves by 204529. In 2020, however, the Mangrove Working Group and its supporting regulations were abolished and the mangrove rehabilitation strategy was subsequently managed by BRGM4. This effectively removed the regional governments (subnational working groups) from decisions related to mangrove management and concentrated development of policy at the level of the national government. The new strategy includes a tenfold increase in the annual rehabilitation target (from 11,250 to ~120,000 ha yr−1) with an overall target of 600,000 ha to be achieved within a shorter timeline (2020–2024). Without clear planning and appropriate strategies, these ambitious targets may not be feasible. For example, the annual mangrove rehabilitation area reached between 2017 and 2020 was only 5,318 ha (50% of the target) despite 2.6 million seedlings being planted (Supplementary Table 3). Given the lessons from the previous mangrove rehabilitation and the emerging processes of mangrove governance, it is timely to set an achievable restoration framework with improved planning, evaluation and monitoring.Implication for international environmental agendasA successful mangrove rehabilitation programme can directly contribute to reducing poverty (SDG 1) and maintaining food security and livelihoods (SDG 2), thereby increasing the health and well-being of 74 million coastal people in Indonesia (see Supplementary Table 1 for total population of regions with restoration potential area >5 ha). Additionally, mangrove rehabilitation will directly contribute to other relevant SDGs, such as improving water quality (SDG 6), providing healthy coastal habitats for fish and other marine biodiversity (SDG 14), contributing to emissions reductions and improving coastal resilience from sea level rise (SDG 13) and sustainably managing and protecting terrestrial ecosystems (SDG 15). Mangrove rehabilitation contributions to SDG 1 and 2 are particularly relevant as the current rehabilitation programme is delivered as cash-for-works activities under the National Economic Recovery strategy (PEN) as part of the social welfare payments to alleviate economic impacts of the COVID-19 pandemic17. With the current annual mangrove rehabilitation budget of US$0.10 billion17, further implementation of scalable community-based mangrove restoration with technical support from subnational and non-government stakeholders could increase the benefits to local communities, if administered properly. Therefore, the large investments planned for coastal communities via a national mangrove restoration programme will not only contribute to the economy of coastal communities, potentially reducing poverty across 199 regencies but will also help in securing nearly 4% of the national greenhouse gas emissions reduction target from the land sector.Restoring 193,367 ha of mangroves in the next 5 years (2021–2025) may contribute to carbon sequestration of 22 ± 10 MtCO2e by 2030 (see Methods for detailed estimate calculation and assumptions). Moreover, stopping the current annual rates of mangrove loss of 7,436 ha yr−1 between 2021 and 2030 will reduce up to 58 ± 37 MtCO2e or 12% of the national land sector emissions reduction targets. Clearly, climate benefits from mangrove rehabilitation and conservation in Indonesia are substantial if rehabilitation and conservation can be implemented appropriately and large annual rehabilitation targets are achieved. Indonesia has submitted its updated Nationally Determined Contributions (NDCs) to the United Nations Framework Convention on Climate Change, within which integrated management and rehabilitation of mangroves is a component of the actions to enhance the resilience of coastal ecosystems30. Further ecological aquaculture practices such as silvofisheries which are commonly applied in Indonesia31,32 may provide promising potential for climate change mitigation through mangrove biomass enhancement. With the increased potential for international investment to support mangrove rehabilitation in Indonesia, there is an opportunity for Indonesia to take the lead and show the world how mangrove conservation and rehabilitation can contribute to multiple international environmental agendas.In the past three decades, the governance of mangrove conservation and rehabilitation in Indonesia has been highly variable in approach (Extended Data Fig. 1). The current approach is top-down4 which has risks and may be ineffective at achieving landscape-scale increases in mangrove extent, as was demonstrated post-tsunami24,29. This top-down approach set by national-level agencies, which are responsible for achieving rehabilitation targets, has limited involvement (or investment) by subnational governments. While we have identified key factors that determine land available for mangrove rehabilitation, the success of mangrove rehabilitation is not necessarily assured because of the limited involvement of subnational mangrove working groups. A current ‘one size fits all’ strategy of the national government may not be appropriate to achieve successful mangrove rehabilitation and thus more flexible, localized approaches may increase the likelihood of success. More

Patterson, J. et al. Nat. Sustain. 4, 841–850 (2021).Article 

Google Scholar 
Reid, A. J. et al. Biol. Rev. 94, 849–873 (2019).Article 

Google Scholar 
Vollmer, D. & Harrison, I. J. Environ. Res. Lett. 16, 011005 (2021).Article 

Google Scholar 
Zeitoun, M. et al. Glob. Environ. Change 39, 143–154 (2016).Article 

Google Scholar 
Bezerra, M. O. et al. Environ. Manage. 69, 815–834 (2022).Article 

Google Scholar 
Souter, N. J. et al. Water 12, 788 (2020).Article 
CAS 

Google Scholar 
Akhmouch, A., Clavreul, D. & Glas, P. Water Int. 43, 5–12 (2018).Article 

Google Scholar 
Andersson, E. Ambio 51, 1–8 (2022).Article 

Google Scholar 
Huntington, H. P. et al. Nat. Sustain. 4, 672–679 (2021).Article 

Google Scholar 
Soames Job, R. F. Am. J. Public Health 78, 163–167 (1988).Article 
CAS 

Google Scholar 
Poff, N. L. et al. Nat. Clim. Change 6, 25–34 (2016).Article 

Google Scholar 
Diaz-Kope, L. & Miller-Stevens, K. Public Works Management and Policy 20, 29–48 (2015).Article 

Google Scholar 
OECD Financing a Water Secure Future (OECD Publishing, 2022).Cardascia, S. Financing Water Infrastructure and Landscape Approaches in Asia and the Pacific. Background Paper for 5th Roundtable on Financing Water (OECD Publishing, 2019).Schlager, E. & Blomquist, W. Embracing Watershed Politics (University Press of Colorado, 2008).Wehn, U., Collins, K., Anema, K., Basco-Carrera, L. & Lerebours, A. Water Int. 43, 34–59 (2018).Article 

Google Scholar 
Shaad, K., Souter, N. J., Vollmer, D., Regan, H. M. & Bezerra, M. O. Environ. Manage. 69, 752–767 (2022).Article 

Google Scholar  More

Potts, S. G. et al. Global pollinator declines: trends, impacts and drivers. Trends Ecol. Evol. 25, 345–353 (2010).Article 
PubMed 

Google Scholar 
Lautenbach, S., Seppelt, R., Liebscher, J. & Dormann, C. F. Spatial and temporal trends of global pollination benefit. PLoS ONE 7, e35954 (2012).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Ollerton, J., Winfree, R. & Tarrant, S. How many flowering plants are pollinated by animals? Oikos 120, 321–326 (2011).Article 

Google Scholar 
Rodger, J. G. et al. Widespread vulnerability of flowering plant seed production to pollinator declines. Sci. Adv. 7, eabd3524 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Biesmeijer, J. C. et al. Parallel declines in pollinators and insect-pollinated plants in Britain and the Netherlands. Science 313, 351–354 (2006).Article 
CAS 
PubMed 

Google Scholar 
Bennett, J. M. et al. Land use and pollinator dependency drives global patterns of pollen limitation in the Anthropocene. Nat. Commun. 11, 3999 (2020).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Tylianakis, J. M., Didham, R. K., Bascompte, J. & Wardle, D. A. Global change and species interactions in terrestrial ecosystems. Ecol. Lett. 11, 1351–1363 (2008).Article 
PubMed 

Google Scholar 
Hegland, S. J., Nielsen, A., Lázaro, A., Bjerknes, A.-L. & Totland, Ø. How does climate warming affect plant–pollinator interactions? Ecol. Lett. 12, 184–195 (2009).Article 
PubMed 

Google Scholar 
Thébault, E. & Fontaine, C. Stability of ecological communities and the architecture of mutualistic and trophic networks. Science 329, 853–856 (2010).Article 
PubMed 

Google Scholar 
Lever, J. J., van Nes, E. H., Scheffer, M. & Bascompte, J. The sudden collapse of pollinator communities. Ecol. Lett. 17, 350–359 (2014).Article 
PubMed 

Google Scholar 
Valdovinos, F. S. et al. Species traits and network structure predict the success and impacts of pollinator invasions. Nat. Commun. 9, 2153 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Waser, N. M., Chittka, L., Price, M. V., Williams, N. M. & Ollerton, J. Generalization in pollination systems, and why it matters. Ecology 77, 1043–1060 (1996).Article 

Google Scholar 
Brosi, B. J. Pollinator specialization: from the individual to the community. New Phytol. 210, 1190–1194 (2016).Article 
PubMed 

Google Scholar 
Elmqvist, T. et al. Response diversity, ecosystem change, and resilience. Front. Ecol. Environ. 1, 488–494 (2003).Article 

Google Scholar 
Waser, N. M. & Ollerton, J. Plant–Pollinator Interactions: From Specialization to Generalization (Univ. of Chicago Press, 2006).Ashman, T.-L., Arceo-Gómez, G., Bennett, J. M. & Knight, T. M. Is heterospecific pollen receipt the missing link in understanding pollen limitation of plant reproduction? Am. J. Bot. 107, 845–847 (2020).Article 
PubMed 

Google Scholar 
Garibaldi, L. A. et al. Trait matching of flower visitors and crops predicts fruit set better than trait diversity. J. Appl. Ecol. 52, 1436–1444 (2015).Article 

Google Scholar 
CaraDonna, P. J. et al. Seeing through the static: the temporal dimension of plant–animal mutualistic interactions. Ecol. Lett. 24, 149–161 (2021).Article 
PubMed 

Google Scholar 
Burkle, L. A., Marlin, J. C. & Knight, T. M. Plant–pollinator interactions over 120 years: loss of species, co-occurrence, and function. Science 339, 1611–1615 (2013).Article 
CAS 
PubMed 

Google Scholar 
Jacquemin, F. et al. Loss of pollinator specialization revealed by historical opportunistic data: insights from network-based analysis. PLoS ONE 15, e0235890 (2020).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Mathiasson, M. E. & Rehan, S. M. Wild bee declines linked to plant–pollinator network changes and plant species introductions. Insect Conserv. Divers. 13, 595–605 (2020).Article 

Google Scholar 
Bennett, J. M. et al. A review of European studies on pollination networks and pollen limitation, and a case study designed to fill in a gap. AoB Plants 10, ply068 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Doré, M., Fontaine, C. & Thébault, E. Relative effects of anthropogenic pressures, climate, and sampling design on the structure of pollination networks at the global scale. Glob. Change Biol. 27, 1266–1280 (2021).Article 

Google Scholar 
Rader, R. et al. Non-bee insects are important contributors to global crop pollination. Proc. Natl Acad. Sci. USA 113, 146–151 (2016).Article 
CAS 
PubMed 

Google Scholar 
Post, E. et al. Ecological dynamics across the arctic associated with recent climate change. Science 325, 1355–1358 (2009).Article 
CAS 
PubMed 

Google Scholar 
Hung, K.-L. J., Kingston, J. M., Albrecht, M., Holway, D. A. & Kohn, J. R. The worldwide importance of honey bees as pollinators in natural habitats. Proc. R. Soc. B 285, 20172140 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Kearns, C. A. Anthophilous fly distribution across an elevation gradient. Am. Midl. Nat. 127, 172–182 (1992).Article 

Google Scholar 
Kevan, P. G. Insect pollination of high arctic flowers. J. Ecol. 60, 831–847 (1972).Article 

Google Scholar 
Tiusanen, M., Hebert, P. D. N., Schmidt, N. M. & Roslin, T. One fly to rule them all—muscid flies are the key pollinators in the arctic. Proc. Roy. Soc. B 283, 20161271 (2016).Article 

Google Scholar 
Weiner, C., Werner, M., Linsenmair, K. E. & Blüthgen, N. Land use intensity in grasslands: changes in biodiversity, species composition and specialisation in flower visitor networks. Basic Appl. Ecol. 12, 292–299 (2011).Article 

Google Scholar 
Rader, R., Edwards, W., Westcott, D. A., Cunningham, S. A. & Howlett, B. G. Pollen transport differs among bees and flies in a human-modified landscape. Divers. Distrib. 17, 519–529 (2011).Article 

Google Scholar 
Bartley, T. J. et al. Food web rewiring in a changing world. Nat. Ecol. Evol. 3, 345–354 (2019).Article 
PubMed 

Google Scholar 
Ghisbain, G., Gérard, M., Wood, T. J., Hines, H. M. & Michez, D. Expanding insect pollinators in the Anthropocene. Biol. Rev. 96, 2755–2770 (2021).Article 
PubMed 

Google Scholar 
Silén, F. Blombiologiska iakttagelser i Kittilä Lappmark. Medd. Soc. Fauna Flora Fennica 31, 80–99 (1906).
Google Scholar 
Clavel, J., Julliard, R. & Devictor, V. Worldwide decline of specialist species: toward a global functional homogenization? Front. Ecol. Environ. 9, 222–228 (2011).Article 

Google Scholar 
Erhardt, A. Pollination of Dianthus superbus L. Flora 185, 99–106 (1991).Article 

Google Scholar 
Witt, T., Jürgens, A., Geyer, R. & Gottsberger, G. Nectar dynamics and sugar composition in flowers of Silene and Saponaria species (Caryophyllaceae). Plant Biol. 1, 334–345 (1999).Article 
CAS 

Google Scholar 
Morales, C. L. & Traveset, A. Interspecific pollen transfer: magnitude, prevalence and consequences for plant fitness. Crit. Rev. Plant Sci. 27, 221–238 (2008).Article 
CAS 

Google Scholar 
Ashman, T.-L. & Arceo-Gómez, G. Toward a predictive understanding of the fitness costs of heterospecific pollen receipt and its importance in co-flowering communities. Am. J. Bot. 100, 1061–1070 (2013).Article 
PubMed 

Google Scholar 
Orford, K. A., Vaughan, I. P. & Memmott, J. The forgotten flies: the importance of non-syrphid Diptera as pollinators. Proc. R. Soc. B 282, 20142934 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Stavert, J. R. et al. Hairiness: the missing link between pollinators and pollination. PeerJ 4, e2779 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Doyle, T. et al. Pollination by hoverflies in the Anthropocene. Proc. R. Soc. B 287, 20200508 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Albrecht, M., Schmid, B., Hautier, Y. & Müller, C. B. Diverse pollinator communities enhance plant reproductive success. Proc. R. Soc. B. 279, 4845–4852 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
Fründ, J., Dormann, C. F., Holzschuh, A. & Tscharntke, T. Bee diversity effects on pollination depend on functional complementarity and niche shifts. Ecology 94, 2042–2054 (2013).Article 
PubMed 

Google Scholar 
Magrach, A., Molina, F. P. & Bartomeus, I. Niche complementarity among pollinators increases community-level plant reproductive success. Peer Commun. J. 1, e1 (2021).Article 

Google Scholar 
Giménez-Benavides, L., Dötterl, S., Jürgens, A., Escudero, A. & Iriondo, J. M. Generalist diurnal pollination provides greater fitness in a plant with nocturnal pollination syndrome: assessing the effects of a Silene–Hadena interaction. Oikos 116, 1461–1472 (2007).
Google Scholar 
Vázquez, D. P., Blüthgen, N., Cagnolo, L. & Chacoff, N. P. Uniting pattern and process in plant–animal mutualistic networks: a review. Ann. Bot. 103, 1445–1457 (2009).Article 
PubMed 
PubMed Central 

Google Scholar 
Vizentin-Bugoni, J., Debastiani, V. J., Bastazini, V. A. G., Maruyama, P. K. & Sperry, J. H. Including rewiring in the estimation of the robustness of mutualistic networks. Methods Ecol. Evol. 11, 106–116 (2020).Article 

Google Scholar 
Brosi, B. J. & Briggs, H. M. Single pollinator species losses reduce floral fidelity and plant reproductive function. Proc. Natl Acad. Sci. USA 110, 13044–13048 (2013).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Pekkarinen, A. & Teräs, I. Zoogeography of Bombus and Psithyrus in northwestern Europe (Hymenoptera, Apidae). Ann. Zool. Fennici 30, 187–208 (1993).
Google Scholar 
Arbetman, M. P., Gleiser, G., Morales, C. L., Williams, P. & Aizen, M. A. Global decline of bumblebees is phylogenetically structured and inversely related to species range size and pathogen incidence. Proc. R. Soc. B 284, 20170204 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Kerr, J. T. et al. Climate change impacts on bumblebees converge across continents. Science 349, 177–180 (2015).Article 
CAS 
PubMed 

Google Scholar 
Arceo-Gómez, G., Barker, D., Stanley, A., Watson, T. & Daniels, J. Plant–pollinator network structural properties differentially affect pollen transfer dynamics and pollination success. Oecologia 192, 1037–1045 (2020).Article 
PubMed 

Google Scholar 
de Santiago-Hernández, M. H. et al. The role of pollination effectiveness on the attributes of interaction networks: from floral visitation to plant fitness. Ecology 100, e02803 (2019).Article 
PubMed 

Google Scholar 
Koch, V., Zoller, L., Bennett, J. M. & Knight, T. M. Pollinator dependence but no pollen limitation for eight plants occurring north of the Arctic Circle. Ecol. Evol. 10, 13664–13672 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Loboda, S., Savage, J., Buddle, C. M., Schmidt, N. M. & Høye, T. T. Declining diversity and abundance of High Arctic fly assemblages over two decades of rapid climate warming. Ecography 41, 265–277 (2018).Article 

Google Scholar 
Høye, T. T., Post, E., Schmidt, N. M., Trøjelsgaard, K. & Forchhammer, M. C. Shorter flowering seasons and declining abundance of flower visitors in a warmer Arctic. Nat. Clim. Change 3, 759–763 (2013).Article 

Google Scholar 
Soroye, P., Newbold, T. & Kerr, J. Climate change contributes to widespread declines among bumble bees across continents. Science 367, 685–688 (2020).Article 
CAS 
PubMed 

Google Scholar 
Zattara, E. E. & Aizen, M. A. Worldwide occurrence records suggest a global decline in bee species richness. One Earth 4, 114–123 (2021).Article 

Google Scholar 
Bartomeus, I., Stavert, J. R., Ward, D. & Aguado, O. Historical collections as a tool for assessing the global pollination crisis. Philos. Trans. R. Soc. B 374, 20170389 (2019).Article 

Google Scholar 
Rakosy, D., Ashman, T.-L., Zoller, L., Stanley, A. & Knight, T. M. Integration of historic collections can shed light on patterns of change in plant–pollinator interactions and pollination service. Func. Ecol. https://doi.org/10.1111/1365-2435.14211 (2022).Hyne, C. J. C. W. Through Arctic Lapland (A. and C. Black, 1898).Knuth, P. Handbuch der Blütenbiologie, unter Zugrundelegung von Herman Müllers Werk: ‘Die Befruchtung der Blumen durch Insekten’ (W. Engelmann, 1898).Zoller, L. & Knight, T. M. Historical records of plant-insect interactions in subarctic Finland.BMC Res. Notes 15, 317 (2022).Article 
PubMed 
PubMed Central 

Google Scholar 
Zoller, L. & Knight, T. M. Historical records of plant–insect interactions in subarctic Finland. figshare https://doi.org/10.6084/m9.figshare.c.5828663.v4 (2022).Zoller, L., Bennett, J. M. & Knight, T. M. Diel-scale temporal dynamics in the abundance and composition of pollinators in the arctic summer. Sci. Rep. 10, 21187 (2020).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2020).Hsieh, T. C., Ma, K. H. & Chao, A. iNEXT: an R package for rarefaction and extrapolation of species diversity (Hill numbers). Methods Ecol. Evol. 7, 1451–1456 (2016).Article 

Google Scholar 
Klotz, S., Kühn, I. & Durka, W. Biolflor Database (UFZ—Centre for Environmental Research Leipzig-Halle, 2002); https://www.ufz.de/biolflor/index.jspOksanen, J. et al. vegan: Community ecology package. R version 2.5.7 (2020).Chao, A., Chazdon, R. L., Colwell, R. K. & Shen, T.-J. Abundance-based similarity indices and their estimation when there are unseen species in samples. Biometrics 62, 361–371 (2006).Article 
PubMed 

Google Scholar 
Dormann, C. F. et al. bipartite: Visualising bipartite networks and calculating some (ecological) indices. R version 2.16 (2021).Blüthgen, N., Menzel, F. & Blüthgen, N. Measuring specialization in species interaction networks. BMC Ecol. 6, 9 (2006).Article 
PubMed 
PubMed Central 

Google Scholar 
Stefan, V. & Knight, T. M. bootstrapnet: Bootstrap network metrics. R version 1.0.0 https://valentinitnelav.github.io/bootstrapnet/ (2021).Poisot, T., Canard, E., Mouillot, D., Mouquet, N. & Gravel, D. The dissimilarity of species interaction networks. Ecol. Lett. 15, 1353–1361 (2012).Article 
PubMed 

Google Scholar 
Poisot, T. Dissimilarity of species interaction networks: quantifying the effect of turnover and rewiring. Peer Community Journal 2, e35 (2022).Article 

Google Scholar 
Dormann, C. F. How to be a specialist? Quantifying specialisation in pollination networks. Netw. Biol. 1, 1 (2011).
Google Scholar  More

In this section, we introduce flapping flight dynamics and describe the bird model used in our computational framework. Furthermore, we describe how such a dynamical model is used in order to identify steady and level flapping flight regimes, study their stability, and assess their energetic performance.Equations of motion modelling flapping flightFlight dynamics is restricted to the longitudinal plane and thus the bird main body is captured as a rigid-body with three degrees of freedom, i.e. two in translation and one in rotation. This model preserves symmetry with respect to this plane, without any lateral force and moment. The aerodynamic model of the wing relies on the theory of quasi-steady lifting line23. Additionally, the present work does not account for the inertial forces due to the acceleration of the wing, and thus also neglecting the so-called inertial power. This inertial power was shown to be negligible in fast forward flight conditions, in comparison to the other contributions to actuation power24, and is thus systematically neglected in similar work10,11,25 since wing inertia is neglected.The body is thus modelled with a mass (m_b) and a rotational inertia (I_{yy}) about its center of mass. The equations of motion are expressed in the body frame (G(x’, y’, z’)) with unit vectors ((hat{textbf{e}}_{x’}, hat{textbf{e}}_{y’}, hat{textbf{e}}_{z’})), and an origin located at the center of mass, as pictured in Fig. 1a. The state space vector is thus$$begin{aligned} textbf{x} = {u, w, q, theta } end{aligned}$$where u and w are the body velocities along the (x’-) and (z’-)axis and (theta) and q are the pitch angle and its time derivative about the (y’-)axis, respectively. Consequently, the equations of motion read11,13,26$$begin{aligned} begin{aligned} dot{u}&= -qw – gsin theta + frac{1}{m_b}big ( {F_{x’}(textbf{x}(t), t)} \&quad + {F_{x’, t}(textbf{x}(t), t)} + D (textbf{x}(t), t) big )\ dot{w}&= qu + gcos theta +frac{1}{m_b} big ( F_{z’}(textbf{x}(t), t) + F_{z’, t}(textbf{x}(t), t) big ) \ dot{q}&=frac{1}{I_{yy}} big ( M_{y’}(textbf{x}(t), t) + M_{y’, t}(textbf{x}(t), t) big )\ dot{theta }&= q end{aligned} end{aligned}$$
(1)
Figure 1(a) Bird model for describing the flight dynamics in the longitudinal plane. The state variables are expressed with respect to the moving body-frame located at the flier’s center of mass (G(x’,z’)). These state variables are the component of forward flight velocity, u, the velocity component of local vertical velocity, w, the orientation of this body-centered moving frame with respect to the fixed frame, (theta) and its angular velocity, q. A second frame (O(x’_{w}, z’_{w})) is used to compute the position of the wing, relative to the body. The wings (dark gray) and the tail (red) are the surfaces of application of aerodynamic forces. (b) Top view of the bird model. The left wing emphasizes a cartoon model of the skeleton. The shoulder joint s connects the wing to the body via three rotational degrees of freedom (RDoF), the elbow joint e connects the arm with the forearm via one RDoF and the wrist joint w connects the forearm to the hand via two RDoF. Each feather is attached to a bone via two additional RDoF, except the most distal one ”1” which is rigidly aligned with the hand. The right wing further emphasizes the lifting line (red) which is computed as a function of the wing morphing. The aerodynamic forces generated on the wing are computed on the discretized elements (P_{i}). The tail is modelled as a triangular shape with fixed chord (c_{t}) and maximum width (b_{t}) that can be morphed as a function of its opening angle (beta). (c) Wing element i between two wing profiles, identifying a plane (Sigma) containing the lifting line (red). (d) Cross section of the wing element, containing the chord point (mathbf {P_i}) where the velocities are computed (Color figure online).Full size imageThe forcing terms in Eq. (1) are the aerodynamic forces and moments applied to the wing (namely (F_{x’}), (F_{z’}), and (M_{y’}) ) and to the tail ((F_{x’, t}), (F_{z’, t}), and (M_{y’, t})). The whole drag is captured by an extra force D that sums contributions due to the body (D_{b}), the skin friction of the wing (wing profile) (D_{p,w}), and the skin friction of the tail (tail profile) (D_{p,t}). These terms are described in detail in the next sections. Importantly, we accounted for the drag acting purely along (x’) direction, after proving that the projection of the drag forces along (z’)-axis is between two and three orders of magnitude smaller with respect to the aerodynamic forces produced by two other main lifting surfaces. This assumptions holds for the fast forward flight regime that are subject of our study, but such components of drag along (z’) axis should be accounted for other flight situations.Wing modelThe bird has two wings. Each wing is a rigid poly-articulated body, comprising the bird arm, forearm and hand, as pictured in Fig. 1b. Each segment is actuated by a joint to induce wing morphing. We refer to13,15 for a complete description of this wing kinematic model.Each joint is kinematically driven to follow a sinusoidal trajectory specified as:$$begin{aligned} q_{i}(t) = q_{0,i}(t) + A_{i} sin (omega t + phi _{0,i}) end{aligned}$$
(2)
with (omega = 2 pi f) and f being the flapping frequency which is identical for each joint, (q_{0,i}) being the mean angle over a period (or offset), (A_i) the amplitude, and (phi _{0,i}) the relative phase of joint i. A complete wingbeat cycle is therefore described through a set of 19 kinematic parameters, including the frequency f.We assume that the wing trajectory is rigidly constrained, and therefore we do not need to explicitly solve the wing dynamics. Under this assumption, the motion generation does not require the computation of joint torques. The model further embeds seven feathers of length (l_{ki}) in each wing. The feathers in the model have to be considered a representative sample of the real wing feathers. They thus have a limited biological relevance; their number is chosen so as to interpolate the planform satisfactorily and to smoothly capture the morphing generated by the bone movements. These feathers are attached to their respective wing bones via two rotational degrees of freedom allowing them to pitch and spread in the spanwise direction. These two degrees of freedom are again kinematically driven by relationships that depend on the angle between the wing segments13. This makes the feathers spreading and folding smoothly through the wingbeat cycle. In sum, the kinematic model of the wing yields the position of its bones and feathers at every time step. This provides a certain wing morphing from which the wing envelope (leading edge and trailing edge) can be computed (see Fig. 1b). From the wing envelope, the aerodynamic chord and the lifting line are computed. The lifting line is the line passing through the quarter of chord, which is itself defined as the segment connecting the leading edge to the trailing edge and orthogonal to the lifting line (Fig. 1b). This extraction algorithm is explained in detail in15.In order to calculate the aerodynamic forces, the angle of attack of the wing profile has to be evaluated. Each wing element defines a plane containing the lifting line and the aerodynamic chord as pictured in Fig. 1c. The orientation of the plane is identified by the orthogonal unit vectors ((hat{textbf{e}}_n, hat{textbf{e}}_t, hat{textbf{e}}_b)), where (hat{textbf{e}}_n) is the vector perpendicular to the plane and (hat{textbf{e}}_t) is the tangent to the lifting line. To compute the effective angle of attack, the velocity perceived by the wing profile is computed as the sum of the velocities due to the body and wing motion, and the velocity induced by the wake. The first contribution, (textbf{U}), accounts for$$begin{aligned} textbf{U} = textbf{U}_{infty } – textbf{v}_{kin} – textbf{v}_{q}end{aligned}$$where (textbf{U}_{infty } = u hat{textbf{e}}_{x’} + w hat{textbf{e}}_{z’}) is the actual flight velocity, (textbf{v}_{kin}) is the relative velocity of the wing due to its motion, and (textbf{v}_{q}) is the component induced by the angular velocity of the body q and calculated as$$begin{aligned} textbf{v}_{q} = qhat{textbf{e}}_{y’} wedge (textbf{P}_{i} – textbf{G})end{aligned}$$This velocity vector (textbf{U}) defines the angle (alpha), as pictured in Fig. 1d.The second contribution is due to the induced velocity field by the wake, i.e. the downwash velocity (w_{d}), and acting along the normal unit vector (-w_{d}hat{textbf{e}}_n). The resulting effective angle of attack, (alpha _{r}), is thus$$begin{aligned} alpha _{r} = alpha – frac{w_{d}}{|textbf{U}|}end{aligned}$$The downwash velocity (w_d) is computed according to the Biot-Savart law23, assuming the wake being shed backwards in the form of straight and infinitely long vortex filaments at each time step of the simulation13,15. This quasi-steady approximation is justified a posteriori by ensuring that our reduced frequency, inversely proportional to the unknown airspeed, never exceeds the value of 0.2, below which the effects of time-dependent wake shapes on wing circulation are negligible (e.g. see discussion in27). Once the downwash is evaluated, it is possible to evaluate the circulation, and consequently the aerodynamic force and moment acting at the element (P_i), i.e. (F_{x’, i}(textbf{x}(t), t), F_{z’, i}(textbf{x}(t), t), M_{y’, i}(textbf{x}(t), t)), as explained in detail in13. We use the thin airfoil theory for the estimation of the lift coefficient, with a slope of (2pi) that saturates at an effective angle of attack (alpha _{r}) of (pm 15^{circ }).Drag production by body and wingThe main body and the wings induce drag that should be accounted for in a model aiming at characterizing energetic performance. Body-induced drag is named parasitic because the body itself does not contribute to lift generation, and only induces skin friction and pressure drag around its envelope28. The total body drag is$$begin{aligned} D_{b} = frac{1}{2}rho C_{d, b} S_{b}|textbf{U}_{infty }|^{2} end{aligned}$$
(3)
where (rho) is the air density. We implemented the model described by Maybury28 to compute the body drag coefficient (C_{d, b}). This depends on the morphology of the bird and the Reynolds number Re according to$$begin{aligned} C_{d,b} = 66.6m_{b}^{-0.511}FR_{t}^{0.9015}S_{b}^{1.063}Re^{-0.197} end{aligned}$$
(4)
with (S_{b}) and (FR_{t}) are respectively the frontal area of the body and the fitness ratio of the bird, and both of them can be estimated from other allometric formulas i.e.28,29.$$begin{aligned} S_{b}= & {} 0.00813m_{b}^{2/3} end{aligned}$$
(5)
$$begin{aligned} FR_{t}= & {} 6.0799m_{b}^{0.1523} end{aligned}$$
(6)
The Reynolds number (Re = rho |textbf{U}_{infty }| overline{c} / mu) is calculated with the reference length of the mean aerodynamic chord (overline{c}), with (mu) being the dynamic viscosity. This model is found to be suitable for Reynolds number in the range (10^{4}-10^{5})28. Another source of drag is the profile drag due to friction between the air and the feathers on the wings. It is the sum of the profile drag at each section along the wingspan, i.e.$$begin{aligned} D_{p,w} = frac{1}{2} rho C_{d, pro} sum _{i=1}^{n} c_{i}|textbf{U}_{r,i}|^{2} ds_{i} end{aligned}$$
(7)
with (c_{i}) the chord length, (ds_{i}) the length of the lifting line element along the wingspan, and (textbf{U}_{r,i}) the velocity at the wing section i accounting for the body velocity, the kinematics velocity of the wing and the downwash velocity (Fig. 1c,d). We used a value of profile drag of (C_{d, pro} = 0.02) and this is assumed to be constant over the wingspan and throughout the flapping cycle30. In reality, due to the wing motion, this value should be gait dependent. However, the aforementioned assumption has been largely used in previous works31,32.Tail modelSince the span of the tail is of the same magnitude as its aerodynamic chord, here the lifting line approach cannot be used23. Therefore, the tail is modelled according to the slender delta wing theory, as a triangular planform33. Its morphology is illustrated in Fig. 1b and characterized by the opening angle (beta) and the chord (c_t). This latter parameter is kept constant, thus the tail span is controlled via (beta) from the trigonometrical relationship$$begin{aligned} b_{t} = 2c_{t}tan frac{beta }{2}end{aligned}$$The main limitation of this framework is the low range of angles of attack ((alpha _{tail} More