Ecology

We used high-throughput qPCR to screen for 58 infective agents in four Atlantic salmon farm cohorts from British Columbia throughout their production cycles. We measured presence and copy number for target genetic sequences, characteristic of specific viral, bacterial, and eukaryotic agents, including several recently discovered viruses30,31, known or suspected to cause disease in salmon. These agents displayed various temporal patterns of prevalence and intensity, with several displaying elevated levels in dead and dying fish.
The data and analyses we have presented provide a unique look into the epidemiology of farmed salmon populations, and wildlife/livestock diseases generally. No past studies have had access to multiple farmed-salmon cohorts, throughout their production cycles, with the capacity to molecularly screen for a large suite of infectious agents. Other work has reported agent data for dead-sampled fish collected in BC as part of Fisheries and Oceans Canada’s farm audit program40, but such analyses lack the time-series nature of the results we have presented. To our knowledge, no other studies have provided such detailed, comprehensive information for infective agents in domestic or wild populations over time. This study, therefore, presents a substantial step toward effectively monitoring shared wildlife/livestock diseases, made possible by cutting-edge technology, as predicted previously22.
While our findings offer specific insight to salmon farmers, aquaculture managers, and those concerned with the disease ecology of sympatric wild salmon, we caution that results remain correlative, and relevant patterns (e.g. apparent mortality signatures) require further investigation. Unfortunately, a lack of regular testing of dead and dying fish (collection was opportunistic, at farms’ discretion) resulted in potential for associated patterns to be obscured. Due to this potential bias in the sampling design, we are unable to draw conclusions related to farm-level mortality rates, but several agents showed patterns of note, including elevated levels in dead and dying fish.
Agent patterns
Perhaps the clearest single-agent pattern—the elevated load of T. maritimum in dead and dying fish (Fig. 3B)—matches generally accepted patterns in aquaculture. Induced tenacibaculosis can be responsible for substantial on-farm mortality worldwide41, and mouthrot resulting from T. maritimum in the east Pacific causes substantial losses42. In our study, mouthrot was noted during veterinarians’ sample processing for cohorts one, three, and four in the months after ocean entry. We note that elevated levels in dead and dying fish could represent the bacterium’s acknowledged role as an opportunistic pathogen41, rather than a direct cause of mortality. We also note the positive correlations between T. maritimum load and that of a number of different agents (Fig. 5), consistent with the view that T. maritimum may facilitate co-infections43. Given its high overall prevalence in fish (Table 3), secondary factors—such as co-infections—might exacerbate infection with T. maritimum, playing a role in mortality.
K. thyrsites intensity was elevated in dead and dying fish for cohorts three and four, around the time they were transferred to their final marine locations (Fig. S8). In both cases, the cohorts finished their production cycles in farms on the central east coast of Vancouver Island (Fig. 1), a region in which the risk of K. thyrsites infection is acknowledged to be high44. This myxozoan parasite is economically important due to post-mortem myoliquefaction seen in infected fish, but it is not considered a pathogen45, and it is unclear why higher gene-copy levels would be observed in dead/dying fish. K. thyrsites may merely replicate faster in stressed fish (in this case due to transport). Our surveillance of pathogens did not include skeletal muscle tissue, where K. thyrsites spores develop, but it did include heart, which can be infected by the parasite46. We note that K. thyrsites was correlated with PRV, with both agents known to infect muscle tissue (although red blood cells are the primary infective tissue for PRV). Follow-up histopathological investigation may provide some insight into K. thyrsites distribution and any associations with pathology or patterns of co-infection.
PRV, which is the causative agent of Heart and Skeletal Muscle Inflammation (HSMI)47 and has recently generated controversy in BC28,29,48, shows several patterns of note. PRV prevalence increased to near ubiquity over time (Fig. 2D), concurrent with an increase, peak, slight decline, and then stabilisation of intensity (Fig. 3B). Although our perspective is limited to sampled fish, with a noted potential for bias, the observed PRV trends were consistent across all four cohorts, and the intensity patterns are consistent with previously reported dissemination, peak replication, and long-term persistence phases of the virus within hosts29,48. Past findings suggest that PRV may induce an antiviral response in hosts that can protect them against certain co-infections49,50. Perhaps counter to the generality of this claim, PRV and ASCV exhibited the strongest load correlation out of any we observed across our data set (Fig. 4). ASCV was originally isolated from salmon with HSMI, and was initially thought to play a role in the disease51. Other work has found no relationship between ASCV and PRV infections or HSMI52. In the case of a related baitfish calicivirus, however, there is evidence that viral co-infection is linked to disease manifestation53, so further work is needed to tease these relationships apart. In general, dead and dying Atlantic salmon in our study did not show elevated prevalence or intensity of PRV, except shortly after ocean entry in cohort one (Figs. 2D, 3B). This mortality signature corresponds to the onset of lesions diagnostic of HSMI in cohort one, which subsequently spread to affect the majority of that farm population for most of a year29.
The gill chlamydia bacterium, C. Syngnamydia salmonis, showed a consistent trend towards elevated prevalence in dead and dying salmon (Fig. 2B). Observed intensity was low, however, often averaging approximately a single copy (Fig. S15). Sequencing has validated past detections of this agent on the Fluidigm BioMark™, and has also revealed SNP diversity within the primer-binding region, resulting in potential underdetection (Miller et al. unpublished). Given a putative mortality signature and the lack of prior epidemiological study of this recently discovered agent54,55, we would suggest further work on C. Syngnamydia salmonis.
Ephemeral mortality signatures appeared for several agents. F. psychrophilum was clearly elevated in dead and dying fish in-hatchery, although we only had access to two hatchery cohorts and cannot draw general conclusions. Intensity of both the ASCV and CTV-2 were elevated in sampled dead and dying fish from at least one cohort shortly after ocean entry (Figs. 3A, S2). Both viruses were also present in-hatchery. The previously reported Cutthroat trout virus appears to be apathogenic56 in trout, and has been detected in Atlantic salmon57. Little is known about the novel variant for which we screened, although in situ hybridisation has revealed that infection can be systemic and extensive in the brain (Mordecai et al. 2020). As for the ASCV, associated pathology was found to be likely due to PRV contamination51. Extremely limited information about these two viruses, paired with our findings, warrants further investigation (e.g. with histopathology and in situ hybridisation) to determine if either virus is linked with pathology. As both these viruses were detected in Chinook salmon (Mordecai et al. 2020), and considering their high prevalence in Atlantic salmon farms, the potential risk they pose to wild Pacific salmon populations should be a priority for future research.
Infectious agent levels overall, as measured by relative infectious burden, showed a clear trend in smolts coming out of freshwater hatcheries for cohorts three and four. There, infectious burden was much higher (in one case 10 000 times higher) in dead and dying fish than in live-sampled fish. While the effect dissipated once fish entered the marine environment, it is clear that hatchery fish are dying with—or of—elevated levels of infection. The patterns we observed likely reflect the transition from freshwater to saltwater, with a coincident shift in infective-agent communities58. Smoltification has also been associated with immune depression59, and elevated infectious burden around the time of ocean entry may reflect this. Where we had dead/dying hatchery samples, however, infectious burden was elevated weeks before ocean entry, hinting at the potential for problems in-hatchery.
Across all agents, we observed apparent coinfection signals that clearly differed from random chance (Figs. 5, S21). We point out, however, that due to the longitudinal nature of our study and cursory investigation of agent correlations, the correlation results almost certainly indicate shared temporal trends that may or may not indicate underlying interactions. For example, Candidatus Branchiomonas cysticola and Flavobacterium psychrophilum were positively correlated with each other, but negatively correlated with a number of other agents. This could be due to both agents being common in-hatchery but not in the marine environment (Figs. S3, S6), counter to many other agents’ patterns. Alternatively, the correlation could be due to a biological relationship, perhaps in relation to gill health. We intend to follow up on a number of correlational patterns.
Agent idiosyncrasies
Several agents showed unexpected patterns, or patterns that may be connected to their particular biology.
The putative Narna-like virus, a recently discovered agent31, showed elevated prevalence in dead and dying fish (Fig. S9). This pattern was mainly due to over-representation in dead-sampled fish, as we detected the agent in 13.2% of dead fish, 1.6% of moribund fish, and 0.4% of live-sampled fish in saltwater. Given that Narnaviridae, of which the putative Narna-like virus is a member, is thought mainly to infect fungi60, this virus may be associated with a fungal decomposer. This is speculative, however, and recent genomic evidence from across taxa suggests that the Narnaviridae may be much more widespread than previously thought61.
Counter to the common trend, P. pseudobranchicola tended to be less common in dead and dying fish than in live-sampled fish (Fig. 2C), with dead fish, in particular, tending to exhibit the lowest levels (results not shown). P. pseudobranchicola primarily infects the pseudobranch62, a structure near the gills involved in oxygenating blood in the eye. Infection also occurs in tissue collected for this study, especially gill63, and we speculate that loads in dead fish could be reduced due to myxospore release or degradation of delicate gill tissue after host death. Given that we did not sample the pseudobranch, it is likely that our data underestimates the load of this organism.
The sampling environments (freshwater or marine) of several detections were unexpected. In particular, we detected K. thyrsites and T. maritimum (Fig. 3C) in freshwater hatcheries, although these agents are considered marine species64,65. It is possible that these hatcheries introduced saltwater in the weeks before ocean transfer, to prepare smolts for release. We also detected F. psychrophilum, considered a freshwater bacterium66, in marine net pens (Fig. 2A). The bacterium is known to survive in brackish water67, however, and this is not the first time it has been detected in a marine setting40,68.
Broader connections
Not all infective agents cause disease, and even agents that do can be present long before—or long after—clinical symptoms. Our work presents only a piece of the puzzle in what is a multifaceted, complex scenario of shared wildlife/livestock disease in salmon aquaculture. The controversy surrounding PRV in BC, as an example, illustrates this complexity. While conventional lab challenges using PRV from BC sources have failed to reproduce in BC fish the extent of HSMI lesions observed on Norwegian farms48,69, work related to our study has been able to identify and shed light on HSMI, and related jaundice/anemia in Chinook salmon, in BC salmon farms28,29. While we saw a putative mortality signature in one cohort during this study, the normal course of PRV infection was not always associated with mortality (e.g. Figs. 2D, 3B). More work will be required to elucidate the nuances of PRV infection, factors that induce associated disease, and possible resultant mortality. A fruitful place to start would be to carry out sampling and diagnostics of dead and dying fish in farms and pens experiencing elevated mortality.
Although we have shown putative mortality signatures for several infective agents in farmed Atlantic salmon, these are not necessarily the agents that pose the greatest risk to wild salmon. For one thing, a given agent need not produce the same effects in different species28,70. For another, contact between populations may not coincide with infection maxima. Depending on when farm smolts enter the marine environment, for example, PRV could be at low prevalence in the spring, when a number of wild Pacific salmon species migrate as juveniles15. Other times of year would be more relevant for interactions with other wild species, and there is much scope for transfer between farmed and wild environments. In addressing shared wildlife/livestock disease, we need to consider both wildlife and livestock as populations that serve as potential reservoirs of disease agents, and are susceptible to outbreaks71. In this context, surveillance and monitoring are essential facets of disease management23, providing raw material to develop understanding of disease and build effective management strategies. Parallel work is monitoring wild populations for the same agents we have investigated here, with the prospect of cross-referencing patterns and impacts72,73,74.
The ubiquity of infectious agents on the farms leads naturally to discussion of potential control strategies, which present a variety of challenges in aquaculture. Vaccination has proven successful at times, but the salmon aquaculture industry has a somewhat chequered history with uptake, since vaccination can affect host growth, and thus the bottom line13. In addition, vaccines have only been developed for a handful of agents. Reducing translocations can be an effective control strategy on land20,78, but transmissive properties of the marine environment and highly mobile marine carrier hosts pose challenges to isolating host populations geographically (Krkosek 2015). Our findings provide circumstantial evidence that some agents (e.g. K. thyrsites) respond to translocations. The fact that two of our four focal cohorts moved substantial distances throughout their respective marine production may be cause for concern, considering the infective-agent populations we have shown those cohorts to have harboured. In general, aquaculture-associated disease and related management decisions have a history of generating political controversy75. Infective-agent monitoring and analyses are critical for designing, implementing, and evaluating effective disease-control measures, and for bridging divides in debate surrounding aquaculture.
With respect to the aquaculture industry, the tools we employed in this study may prove useful for disease management and fish health. We have shown that for many agents, patterns of infection in dead and dying fish mirror those in live fish. By integrating high-throughput infectious-agent screening with existing monitoring of dead fish, farm vets and managers could access a wealth of otherwise unavailable or costly information. Combining such results with strategic sacrificial sampling of live fish during mortality peaks could allow additional insight into which agents may be driving mortality. Protecting the ‘herd’ (and its wild neighbours) may justify such mortal sampling. Furthermore, in other ongoing work, we have seen that much infective-agent information is accessible via nonlethal gill biopsy, which also enables high-throughput screening for gene expression patterns associated with various patterns of stress76 and disease77. Used appropriately, such a combination of tools could be very powerful.
Disease monitoring is never complete, and detection always lags behind pathogen spread78, but new technologies—such as those we employed here – can facilitate efficient, lower-cost surveillance and monitoring. Surveillance for existing pathogens and identification of previously unknown pathogens is part of the integrative approach required to understand and control existing and emerging infectious diseases22. Here, we have further demonstrated the utility of high-throughput, modern genetic techniques for monitoring known infective agents and for generating information about previously under-studied agents26,29,30,31,37,38. Further work will target the risk of transfer between wild and farmed hosts and prioritize threats to salmon, farmed and wild. More

1.
Celli, G. & Maccagnani, B. Honey bees as bioindicators of environmental pollution. Bull. Insectol. 56, 1–3 (2003).
Google Scholar 
2.
Quigley, T. P., Amdam, G. V. & Harwood, G. H. Honey bees as bioindicators of changing global agricultural landscapes. Curr. Opin. Insect Sci. 35, 132–137 (2019).
PubMed  Article  Google Scholar 

3.
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 Biol. Sci. 285, 20172140 (2018).
Article  Google Scholar 

4.
Giannini, T. C., Cordeiro, G. D., Freitas, B. M., Saraiva, A. M. & Imperatriz-Fonseca, V. L. The dependence of crops for pollinators and the economic value of pollination in Brazil. J. Econ. Entomol. 108, 849–857 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

5.
Garibaldi, L. A. et al. Mutually beneficial pollinator diversity and crop yield outcomes in small and large farms. Science 351, 388–391 (2016).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

6.
Miñarro, M., García, D. & Martínez-Sastre, R. Los insectos polinizadores en la agricultura: importancia y gestión de su biodiversidad. Ecosistemas Rev. Científica Ecol. y Medio Ambient. 27, 81–90 (2018).
Google Scholar 

7.
Calderone, N. W. Insect pollinated crops, insect pollinators and US agriculture: trend analysis of aggregate data for the period 1992–2009. PLoS ONE 7, 24–28 (2012).
Article  CAS  Google Scholar 

8.
Kaplan, J. K. Colony collapse disorder: an incomplete puzzle. Agric. Res. Mag. 60, 2489 (2012).
Google Scholar 

9.
VanEngelsdorp, D. et al. Colony collapse disorder: a descriptive study. PLoS ONE 4, 1–17 (2009).
Article  CAS  Google Scholar 

10.
Sanchez-Bayo, F. & Goka, K. Pesticide residues and bees: a risk assessment. PLoS ONE 9, e94482 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

11.
Sindiveg, S. N. da I. de P. para D. V. Mapeamento De Abelhas Participativo (MAP). Relatório 3 anos (2014-2017). Colmeia Viva (2017). 61p. https://www.colmeiaviva.com.br/wp-content/uploads/2019/10/RelatorioMAP.pdf

12.
Wolff, L. F. & Santos, R. S. S. Abelhas melíferas: bioindicadores de qualidade ambiental e de sustentabilidade da agricultura familiar de base ecológica. Embrapa Clima Temperado-Documentos https://www.infoteca.cnptia.embrapa.br/infoteca/bitstream/doc/543990/1/documento244.pdf (2008).

13.
Amaro, P. & Godinho, J. Pesticidas e abelhas. Rev. Ciências Agrárias 35, 53–62 (2012).
Google Scholar 

14.
Tosi, S. & Nieh, J. C. A common neonicotinoid pesticide, thiamethoxam, alters honey bee activity, motor functions, and movement to light. Sci. Rep. 7, 1–13 (2017).
Article  CAS  Google Scholar 

15.
Catae, A. F. et al. MALDI-imaging analyses of honeybee brains exposed to a neonicotinoid insecticide. Pest Manag. Sci. 75, 607–615 (2019).
CAS  PubMed  Article  Google Scholar 

16.
Alves, S. Controle Microbiano de Insetos (FEALQ, Piracicaba, 1998).
Google Scholar 

17.
Hajek, A. E. & Eilenberg, J. Natural Enemies: An Introduction to Biological Control (Cambridge University Press, Cambridge, 2018).
Google Scholar 

18.
Lacey, L. A. et al. Insect pathogens as biological control agents: back to the future. J. Invertebr. Pathol. 132, 1–41 (2015).
CAS  PubMed  Article  Google Scholar 

19.
Melo, A. L. D. A., Soccol, V. T. & Soccol, C. R. Bacillus thuringiensis: mechanism of action, resistance, and new applications: a review. Crit. Rev. Biotechnol. 36, 317–326 (2014).
PubMed  Article  CAS  PubMed Central  Google Scholar 

20.
Eski, A., Demir, İ, Sezen, K. & Demirbağ, Z. A new biopesticide from a local Bacillus thuringiensis var. tenebrionis (Xd3) against alder leaf beetle (Coleoptera: Chrysomelidae). World J. Microbiol. Biotechnol. 33, 35 (2017).
Article  CAS  Google Scholar 

21.
Göktürk, T. Pyrethrum ve Bacillus thuringiensis biyopestisitlerinin Pristiphora abietina (Christ, 1791) (Hymenoptera: Tenthredinidae) üzerindeki etkisi. Artvin Çoruh Üniversitesi Orman Fakültesi Derg. 18, 83–87 (2017).
Article  Google Scholar 

22.
MAPA, M. da A. P. e A. Agrofit. MAPA http://agrofit.agricultura.gov.br/agrofit_cons/principal_agrofit_cons (2019).

23.
Konecka, E., Kaznowski, A., Stachowiak, M. & Maciąg, M. Activity of spore-crystal mixtures of new Bacillus thuringiensis strains against Dendrolimus pini (Lepidoptera: Lasiocampidae) and Spodoptera exigua (Lepidoptera: Noctuidae). Folia For. Pol. Ser. A 60, 91–98 (2018).
Google Scholar 

24.
Gupta, S. & Dikshit, A. K. Biopesticides: an ecofriendly approach for pest control. J. Biopestic. 3, 186–188 (2010).
Google Scholar 

25.
Habid, M. E. M. & Andrade, C. F. S. Bactérias entomopatogênicas. in Controle microbiano de insetos (ed Alves, S. B.) 384–446 (FEALQ, Piracicaba, 1998).

26.
Bravo, A., Gill, S. S. & Soberón, M. Mode of action of Bacillus thuringiensis Cry and Cyt toxins and their potential for insect control. Toxicon 49, 423–435 (2007).
CAS  PubMed  Article  PubMed Central  Google Scholar 

27.
Slessor, K. N., Winston, M. L. & Conte, Y. L. Pheromone communication in the honeybee (Apis mellifera L.). J. Chem. Ecol. 31, 2731–2745 (2005).
CAS  PubMed  Article  PubMed Central  Google Scholar 

28.
Lugtenberg, B. Principles of Plant-Microbe Interactions (Springer, Berlin, 2015).
Google Scholar 

29.
Palma, L. & Berry, C. Understanding the structure and function of Bacillus thuringiensis toxins. Toxicon 109, 1–3 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

30.
Malone, L. A. et al. Effects of ingestion of a Bacillus thuringiensis toxin and a trypsin inhibitor on honey bee flight activity and longevity. Apidologie 32, 57–68 (2001).
CAS  Article  Google Scholar 

31.
Yi, D., Fang, Z. & Yang, L. Effects of Bt cabbage pollen on the honeybee Apis mellifera L. Sci. Rep. 8, 1–6 (2018).
Article  CAS  Google Scholar 

32.
D’Urso, V. et al. Observations on midgut of Apis mellifera workers (Hymenoptera: Apoidea) under controlled acute exposures to a Bacillus thuringiensis-based biopesticide. Apidologie 48, 51–62 (2017).
Article  CAS  Google Scholar 

33.
Libardoni, G. et al. Effect of different Bacillus thuringiensis strains on the longevity of Africanized honey bee. Semin. Agrar. 39, 329–338 (2018).
Article  Google Scholar 

34.
Zhu, Y. C., Wang, Y., Portilla, M., Parys, K. & Li, W. Risk and toxicity assessment of a potential natural insecticide, methyl benzoate, in honey bees (Apis mellifera L.). Insects 10, 1–17 (2019).
Google Scholar 

35.
Hesselbach, H., Seeger, J., Schilcher, F., Ankenbrand, M. & Scheiner, R. Chronic exposure to the pesticide flupyradifurone can lead to premature onset of foraging in honeybees Apis mellifera. J. Appl. Ecol. 57, 609–618 (2020).
CAS  Article  Google Scholar 

36.
Gomes, I. N., Vieira, K. I. C., Gontijo, L. M. & Resende, H. C. Honeybee survival and flight capacity are compromised by insecticides used for controlling melon pests in Brazil. Ecotoxicology 29, 97–107 (2020).
Article  CAS  Google Scholar 

37.
Alquisira-Ramírez, E. V., Paredes-Gonzalez, J. R., Hernández-Velázquez, V. M., Ramírez-Trujillo, J. A. & Peña-Chora, G. In vitro susceptibility of Varroa destructor and Apis mellifera to native strains of Bacillus thuringiensis. Apidologie 45, 707–718 (2014).
Article  CAS  Google Scholar 

38.
Fagúndez, G. A., Blettler, D. C., Krumrick, C. G., Bertos, M. A. & Trujillo, C. G. Do agrochemicals used during soybean flowering affect the visits of Apis mellifera L.?. Span. J. Agric. Res. 14, 7 (2016).
Article  Google Scholar 

39.
Alquisira-Ramírez, E. V. et al. Effects of Bacillus thuringiensis strains virulent to Varroa destructor on larvae and adults of Apis mellifera. Ecotoxicol. Environ. Saf. 142, 69–78 (2017).
PubMed  Article  CAS  Google Scholar 

40.
Horta, A. B., Pannuti, L., Baldin, E. L. L. & Furtado, E. L. Toxinas inseticidas de Bacillus thuringiensis. In Biotecnologia Aplicada à Agro&Indústria (ed. Resende, R. R.) 737–773 (Blucher, Erkrath, 2017). https://doi.org/10.5151/9788521211150-21.
Google Scholar 

41.
Malone, L. A., Burgess, E. P. J. & Stefanovic, D. Effects of a Bacillus thuringiensis toxin, two Bacillus thuringiensis biopesticide formulations, and a soybean trypsin inhibitor on honey bee (Apis mellifera L.) survival and food consumption. Apidologie 30, 465–473 (1999).
CAS  Article  Google Scholar 

42.
Potrich, M. et al. Effect of entomopathogens on Africanized Apis mellifera L. (Hymenoptera: Apidae). Rev. Bras. Entomol. 22, 1–2. https://doi.org/10.1016/j.rbe.2017.12.002 (2018).
ADS  Article  Google Scholar 

43.
Wang, Y. Y. et al. Toxicological, biochemical, and histopathological analyses demonstrating that Cry1C and Cry2A are not toxic to larvae of the honeybee Apis mellifera. J. Agric. Food Chem. 63, 6126–6132 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

44.
Renzi, M. T. et al. Chronic toxicity and physiological changes induced in the honey bee by the exposure to fi pronil and Bacillus thuringiensis spores alone or combined. Ecotoxicol. Environ. Saf. 127, 205–213 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

45.
Hendriksma, H. P. et al. Effect of stacked insecticidal cry proteins from maize pollen on nurse bees (Apis mellifera carnica) and their gut bacteria. PLoS ONE 8, 1–11 (2013).
Article  CAS  Google Scholar 

46.
Jia, H. R. et al. The effects of Bt Cry1Ie toxin on bacterial diversity in the midgut of Apis mellifera ligustica (Hymenoptera: Apidae). Sci. Rep. 6, 1–8 (2016).
Article  CAS  Google Scholar 

47.
Jia, H.-R. et al. No effect of Bt Cry1Ie toxin on bacterial diversity in the midgut of the Chinese honey bees, Apis cerana cerana (Hymenoptera, Apidae). Sci. Rep. 7, 1–10 (2017).
Article  CAS  Google Scholar 

48.
Dai, P. et al. The effect of Bt Cry9Ee toxin on honey bee brood and adults reared in vitro, Apis mellifera (Hymenoptera: Apidae). Ecotoxicol. Environ. Saf. 181, 381–387 (2019).
CAS  PubMed  Article  PubMed Central  Google Scholar 

49.
Baptista, A. P. M., Carvalho, G. A., Carvalho, S. M., Carvalho, C. F. & de Bueno Filho, J. S. S. Toxicidade produtos fitossanitários utilizados em citros para Apis mellifera. Ciência Rural 39, 955–961 (2009).
CAS  Article  Google Scholar 

50.
Tomé, H. V. V., Barbosa, W. F., Martins, G. F. & Guedes, R. N. C. Spinosad in the native stingless bee Melipona quadrifasciata: regrettable non-target toxicity of a bioinsecticide. Chemosphere 124, 103–109 (2015).
ADS  PubMed  Article  CAS  Google Scholar 

51.
Kaplan, E. L. & Meier, P. Non parametric estimation from incomplete observation. J. Am. Stat. Assoc. 53, 457–481 (1958).
MATH  Article  Google Scholar 

52.
Therneau, T. M. A Package for Survival Analysis in R. R package version 3.2-7. https://CRAN.R-project.org/package=survival (2020).

53.
Agresti, A. Categorical Data Analysis (Wiley, New York, 2002).
Google Scholar 

54.
Christensen, R. H. B. Ordinal-Regression Models for Ordinal Data. R package version 2019.12-10. https://CRAN.R-project.org/package=ordina (2019).

55.
Russell, L. R Package ’emmeans’: Estimated Marginal Means, aka Least-Squares Means. https://github.com/rvlenth/emmeans (2020). More

1.
Henriksen, J. B. Case study 7. Key principles in implementing ILO convention no. 169. in Research on Best Practices for the Implementation of the Principles of ILO Convention No. 169 (Programme to Promote ILO Convention No. 169, 2008).
2.
MacKay, F. Addressing Past Wrongs: Indigenous Peoples and Protected Areas: The Right to Restitution of Lands and Resources. (Forest Peoples Programme, 2002).

3.
Hunter-gatherers of the Congo Basin: cultures, histories and biology of African Pygmies. (Transaction Publishers, 2014).

4.
Robillard, M. & Bahuchet, S. Les Pygmées et les autres: terminologie, catégorisation et politique. J. Afr. 82, 15–51 (2012).
Google Scholar 

5.
Joiris, D. V. La chasse, la chance, le chant: aspect du systeme ritual des Baka du Cameroun. (Université Libre de Bruxelles, 1998).

6.
Rupp, S. Interethnic Relations in Southeastern Cameroon: Challenging the” Hunter-Gatherer”–” Farmer” Dichotomy. Afr. Study Monogr. 28, 37–56 (2003).
Google Scholar 

7.
Bahuchet, S. Dans la Forêt D’Afrique Centrale: Les Pygmées Aka Et Baka. (PEETERS-SELAF, 1992).

8.
Yasuoka, H. Long-term foraging expeditions (Molongo) among the Baka Hunter-Gatherers in the Northwestern Congo Basin, with Special Reference to the “Wild Yam Question”. Hum. Ecol. 34, 275–296 (2006).
Article  Google Scholar 

9.
Gallois, S. & Duda, R. Beyond productivity: The socio-cultural role of fishing among the Baka of southeastern Cameroon. Revue d’ethnoécologie [En ligne], (2016).

10.
Yasuoka, H. Snare hunting among Baka hunter-gatherers: Implications for sustainable wildlife management. Afr. Study Monogr. 49, 115–136 (2014).
Google Scholar 

11.
Hayashi, K. Hunting activities in forest camps among the Baka hunter-gatherers of Southeastern Cameroon. Afr. Study Monogr. 29, 73–92 (2008).
Google Scholar 

12.
Duda, R., Gallois, S. & Reyes-Garcia, V. Hunting techniques, wildlife offtake and market integration. A perspective from individual variations among the Baka (Cameroon). Afr. Study Monogr. 38, 97–118 (2017).

13.
Pemunta, N. V. Fortress conservation, wildlife legislation and the Baka Pygmies of southeast Cameroon. GeoJournal 84, 1035–1055 (2019).
Article  Google Scholar 

14.
Dounias, E. & Leclerc, C. Spatial shifts and migration time scales among the Baka Pygmies of Cameroon and the Punan of Borneo. in The social ecology of tropical forest: Migration, populations and frontiers (eds. de Jong, W., Lye, T. P. & Ken-Chi, A.) 147–173 (University Press and Trans Pacific Press, 2006).

15.
Leclerc, C. L’adoption de l’agriculture chez les Pygmées Baka du Cameroun. (Editions Quae, 2012).

16.
Gallois, S. Dynamics of Local ecological knowledge among the Baka children from southeastern Cameroon. (Universitat Autònoma de Barcelona, 2015).

17.
Coquery-Vidrovitch, C. Le Congo au temps des grandes compagnies concessionnaires 1898–1930. (Mouton & Co, 2017).

18.
Pourtier, R. Le Gabon: etat et developpement. vol. Tome 2 (l’Harmattan, 1989).

19.
Jean, S. Les jachères en Afrique tropicale: interprétation technique et foncière. (Musée de l’Homme, 1975).

20.
Bailey, R. C., Bahuchet, S. & Hewlett, B. S. Development in the Central African rainforest: concern for forest peoples. in Conservation of West and Central African Rainforest (ed. Cleaver, K. M.) 202–211 (International Union for Conservation of Nature and Natural Resource, 1992).

21.
Bahuchet, S., McKey, D. & de Garine, I. Wild yams revisited: Is independence from agriculture possible for rain forest hunter-gatherers?. Hum. Ecol. 19, 213–243 (1991).
Article  Google Scholar 

22.
Froment, A. Human biology and the health of African rainforest inhabitants. in Hunter-gatherers of the Congo basin 117–164 (Transaction Publishers, 2014).

23.
Althabe, G. Changements sociaux chez les Pygmées Baka de l’Est-Cameroun. Cahiers d’études africaines 561–592 (1965).

24.
Kitanishi, K. Cultivation by the Baka hunter-gatherers in the tropical rain forest of central Africa. Afr. Study Monogr. 28, 143–157 (2003).
Google Scholar 

25.
Knight, J. Relocated to the roadside: preliminary observations on the forest Peoples of Gabon. Afr. Study Monogr. 28, 81–121 (2003).
Google Scholar 

26.
Yasuoka, H. Fledging agriculturalists? Rethinking the adoption of cultivation by the Baka hunter-gatherers. Afr. Study Monogr. 43, 85–114 (2012).
Google Scholar 

27.
Duda, R. Ethnoecology of hunting in an empty forest: practices, local perceptions and social change among the Baka (Cameroon). (Universitat Autonoma de Barcelona, 2017).

28.
Reyes-García, V. et al. Dietary transitions among three contemporary hunter-gatherers across the tropics. Food Secur. 11, 109–122 (2019).
Article  Google Scholar 

29.
Gallois, S., Heger, T., van Andel, T., Sonké, B. & Henry, A. G. From bush mangoes to bouillon cubes: wild plants and diet among the Baka, forager-horticulturalists from Southeast Cameroon. Econ. Bot. https://doi.org/10.1007/s12231-020-09489-x (2020).
Article  Google Scholar 

30.
Fa, J. E. et al. Differences between pygmy and non-pygmy hunting in Congo Basin Forests. PLoS ONE 11, e0161703 (2016).
Article  Google Scholar 

31.
Dounias, E. & Froment, A. When forest-based hunter-gatherers become sedentary: consequences for diet and health. Unasylva 57, 8 (2006).
Google Scholar 

32.
Hagino, I., Sato, H. & Yamauchi, T. The demographic characteristics and nutritional status for a hunter-gatherer society with social transitions in southeastern Cameroon. Afr. Study Monogr. 45–57 (2014).

33.
Mertens, B. & Lambin, E. F. Spatial modelling of deforestation in southern Cameroon: spatial disaggregation of diverse deforestation processes. Appl. Geogr. 17, 143–162 (1997).
Article  Google Scholar 

34.
Mertens, B. & Lambin, E. F. Land-cover-change trajectories in Southern Cameroon. Ann. Assoc. Am. Geogr. 90, 467–494 (2000).
Article  Google Scholar 

35.
Ndoye, O. & Kaimowitz, D. Macro-economics, markets and the humid forests of Cameroon, 1967–1997. J. Mod. Afr. Stud. 38, 225–253 (2000).
Article  Google Scholar 

36.
Geist, H. J. & Lambin, E. F. Proximate causes and underlying driving forces of tropical deforestation. Bioscience 52, 143 (2002).
Article  Google Scholar 

37.
Messerli, P., Giger, M., Dwyer, M. B., Breu, T. & Eckert, S. The geography of large-scale land acquisitions: Analysing socio-ecological patterns of target contexts in the global South. Appl. Geogr. 53, 449–459 (2014).
Article  Google Scholar 

38.
Seaquist, J. W., Li Johansson, E. & Nicholas, K. A. Architecture of the global land acquisition system: Applying the tools of network science to identify key vulnerabilities. Environ. Res. Lett. 9, 114006 (2014).
ADS  Article  Google Scholar 

39.
Yasuoka, H. et al. Changes in the composition of hunting catches in southeastern Cameroon: a promising approach for collaborative wildlife management between ecologists and local hunters. E&S 20, art25 (2015).

40.
Ohl-Schacherer, J. et al. The Sustainability of subsistence hunting by Matsigenka native communities in Manu National Park Peru. Conserv. Biol. 21, 1174–1185 (2007).
Article  Google Scholar 

41.
Parry, L., Barlow, J. & Peres, C. A. Hunting for sustainability in tropical secondary forests. Conserv. Biol. 23, 1270–1280 (2009).
Article  Google Scholar 

42.
Kümpel, N. F., Milner-Gulland, E. J., Rowcliffe, J. M. & Cowlishaw, G. Impact of gun-hunting on diurnal primates in continental Equatorial Guinea. Int. J. Primatol. 29, 1065–1082 (2008).
Article  Google Scholar 

43.
Coad, L. et al. Social and ecological change over a decade in a village hunting system, central Gabon. Conserv. Biol. 27, 270–280 (2013).
CAS  Article  Google Scholar 

44.
Walters, G., Schleicher, J., Hymas, O. & Coad, L. Evolving hunting practices in Gabon: Lessons for community-based conservation interventions. E&S 20, art31 (2015).

45.
Zerca y Lejos. Soberanía alimentaria y medios de vida (Baka food security). https://zercaylejos.org/proyectos/soberania-alimentaria/ (2020).

46.
Zerca y Lejos. Educación. https://zercaylejos.org/proyectos/educacion/ (2020).

47.
Ávila Martin, E. et al. Wild meat hunting and use by sedentarised Baka Pygmies in southeastern Cameroon. PeerJ 8, e9906 (2020).
Article  Google Scholar 

48.
Ndameau, B. Case Study 7: Cameroon – Boumba Bek – Protected areas and indigenous peoples: The paradox of conservation and survival of the Baka in Moloundou region (south-east Cameroon). in Indigenous Peoples and Protected Areas in Africa: From Principles to Practice (eds. Nelson, J. & Hossack, L.) 215–241 (Forest Peoples Programme, 2001).

49.
Egbe, S. E. et al. The Law, Communities and Wildlife Management in Cameroon. in Cameroon Rural Development Forestry Network. Network Paper (Overseas Development Institute, 2001).

50.
World Weather Online. Djoum Monthly Climate Averages. https://www.worldweatheronline.com/djoum-weather-averages/sud/cm.aspx (2020).

51.
Letouzey, R. Notice phytogéographique du Cameroun au 1:500000. (1985).

52.
Global Forest Watch. An overview of logging in Cameroon. (World Resources Institute (WRI), 2000).

53.
Morgan, D. et al. Impacts of selective logging and associated anthropogenic disturbance on intact forest landscapes and apes of Northern Congo. Front. For. Glob. Change 2, 28 (2019).
Article  Google Scholar 

54.
Ichikawa, M. Problems in the conservation of rainforests in Cameroon. Afr. Study Monogr. 33, 3–20 (2006).
Google Scholar 

55.
Ichikawa, M. Forest Conservation and Indigenous Peoples in The Congo Basin: New Trends toward Reconciliation between Global Issues and Local Interest. in Hunter-gatherers of the Congo Basin: cultures, histories and biology of African Pygmies (ed. Hewlett, B. S.) 321–342 (Transaction Publishers, 2014).

56.
Hattori, S. Current issues facing the forest people in southeastern Cameroon: the dynamics of Baka life and their ethnic relationship with farmers. Afr. Study Monogr. 47, 97–119 (2014).
Google Scholar 

57.
United Nations. Free, Prior and Informed Consent of Indigenous Peoples. https://www.ohchr.org/Documents/Issues/ipeoples/freepriorandinformedconsent.pdf (2013).

58.
Darwin Initiative. Project 24029: Enabling Baka attain food security, improved health and sustain biodiversity. (2020).

59.
Remote Pixel. Remote Pixel—Satellite Imagery Search. https://remotepixel.ca/ (2018).

60.
European Space Agency. Sentinel online. https://sentinel.esa.int/web/sentinel/missions/sentinel-2 (2020).

61.
System for Automated Geoscientific Analyses (SAGA). SAGA – System for Automated Geoscientific Analyses. http://www.saga-gis.org/en/index.html (2020).

62.
Global Forest Watch. Cameroon forest management units (last updated 2018). https://data.globalforestwatch.org/datasets/cameroon-forest-management-units (2020).

63.
Global Forest Watch. Cameroon community forests (last updated 2018). http://data.globalforestwatch.org/datasets/cameroon-community-forests?geometry=6.430%2C5.125%2C17.724%2C8.941 (2020).

64.
Calenge, C. & Fortmann-Roe, S. Package ‘adehabitatHR’. (2020).

65.
Benhamou, S. & Cornélis, D. Incorporating movement behavior and barriers to improve kernel home range space use estimates. J. Wildl. Manag. 74, 1353–1360 (2010).
Article  Google Scholar 

66.
Benhamou, S. & Riotte-Lambert, L. Beyond the Utilization Distribution: Identifying home range areas that are intensively exploited or repeatedly visited. Ecol. Model. 227, 112–116 (2012).
Article  Google Scholar 

67.
Lechowicz, M. J. The sampling characteristics of electivity indices. Oecologia 52, 22–30 (1982).
ADS  Article  Google Scholar 

68.
R Foundation for Statistical Computing. R. (2018).

69.
Jarque, C. M. & Bera, A. K. A test for normality of observations and regression residuals. Int. Stat. Rev. 55, 163 (1987).
MathSciNet  Article  Google Scholar 

70.
Gavrilov, I. & Pusev, R. Package ‘normtest’. (2015).

71.
Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. lmerTest package: Tests in linear mixed effects models. J. Stat. Soft. 82, (2017).

72.
Gockowski, J. et al. The forest margins of Cameroon. in Slash-and-burn Agriculture: The Search for Alternatives (eds. Palm, C. A., Vosti, A. S., Sanchez, A. P. & Polly, J. E.) 305–331 (Columbia University Press, 2005).

73.
Nounamo, L. & Yemefack, M. Shifting cultivation practices and sustainable forestland use in the evergreen forest of Cameroon. in Seminar proceedings ‘Sustainable management of African rain forest’, held inKribi, Cameroon, November 1999. Part II. (eds. Jonkers, W. B. J., Foahom, B. & Schmidt, P.) (The Tropenbos Foundation, 2001).

74.
Chapin, M., Lamb, Z. & Threlkeld, B. Mapping indigenous lands. Annu. Rev. Anthropol. 34, 619–638 (2005).
Article  Google Scholar 

75.
Rainforest Foundation UK. Mapping for rights: Putting communities on the map. Mapping for rights. https://www.mappingforrights.org/index (2020).

76.
Fisher, R., Heckbert, S., León Villalobos, J. M. & Sutton, S. Augmenting physical 3D models with projected information to support environmental knowledge exchange. Appl. Geogr. 112, 102095 (2019).
Article  Google Scholar 

77.
Njounan Tegomo, O., Defo, L. & Usongo, L. Mapping of resource use area by the Baka Pygmies inside and around Boumba-Bek National Park in Southeast Cameroon, with special reference to Baka’s customary rights. Afr. Study Monogr. 43, 45–59 (2012).

78.
Bobo, K. S., Kamgaing, T. O. W., Kamdoum, E. C. & Dzefack, Z. C. B. Bushmeat hunting in southeastern Cameroon: magnitude and impact on duikers (Cephalophus spp.). African Study Monographs Suppl. 51, 119–141 (2015).

79.
Constantino, P. de A. L. Dynamics of hunting territories and prey distribution in Amazonian Indigenous Lands. Appl. Geogr. 56, 222–231 (2015).

80.
Read, J. M. et al. Space, place, and hunting patterns among Indigenous Peoples of the Guyanese Rupununi region. Journal of Latin American Geography 9, 213–243 (2010).
Article  Google Scholar 

81.
Smith, D. A. The spatial patterns of indigenous wildlife use in western Panama: Implications for conservation management. Biol. Cons. 141, 925–937 (2008).
Article  Google Scholar 

82.
Pangau-Adam, M., Noske, R. & Muehlenberg, M. Wildmeat or bushmeat? Subsistence hunting and commercial harvesting in Papua (West New Guinea) Indonesia. Hum. Ecol. 40, 611–621 (2012).
Article  Google Scholar 

83.
Kümpel, N. F. Incentives for sustainable hunting of bushmeat in Río Muni, Equatorial Guinea. (Imperial college London, 2006).

84.
Benítez-López, A. et al. The impact of hunting on tropical mammal and bird populations. Science 356, 180–183 (2017).
ADS  Article  Google Scholar 

85.
Benítez-López, A., Alkemade, R. & Verweij, P. A. The impacts of roads and other infrastructure on mammal and bird populations: A meta-analysis. Biol. Cons. 143, 1307–1316 (2010).
Article  Google Scholar 

86.
World Land Trust. Area converter. https://www.worldlandtrust.org/get-involved/educational-resources/area-converter/ (2020).

87.
Billong Fils, P. E. et al. Ethnobotanical survey of wild edible plants used by Baka People in southeastern Cameroon. J. Ethnobiol. Ethnomed. 16 (2020).

88.
Garnett, S. T. et al. A spatial overview of the global importance of Indigenous lands for conservation. Nat. Sustain. 1, 369–374 (2018).
Article  Google Scholar 

89.
Bruce, T. et al. Faunal Inventory of the Dja Faunal Reserve, Cameroon—2018. 62 (2018). More

1.
Bruno, J. F. & Valdivia, A. Coral reef degradation is not correlated with local human population density. Sci. Rep. 6, 1–8 (2016).
Article  CAS  Google Scholar 
2.
Hughes, T. P. et al. Global warming transforms coral reef assemblages. Nature 556, 492–496 (2018).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

3.
Mollica, N. R. et al. Ocean acidification affects coral growth by reducing skeletal density. Proc. Natl. Acad. Sci. 115, 1754–1759 (2018).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

4.
Mumby, P. J. Stratifying herbivore fisheries by habitat to avoid ecosystem overfishing of coral reefs. Fish Fish. 17, 266–278 (2016).
Article  Google Scholar 

5.
Silbiger, N. J. et al. Nutrient pollution disrupts key ecosystem functions on coral reefs. Proc. R. Soc. B Biol. Sci. 285, 1–9 (2018).
Google Scholar 

6.
Hughes, T. P. Catastrophes, phase shifts, and large-scale degradation of a caribbean coral reef. Science (80-. ). 265, 1547–1551 (1994).

7.
Pandolfi, J. M. et al. Global trajectories of the long-term decline of coral reef ecosystems. Science 301, 955–958 (2003).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

8.
Spalding, M. D. & Brown, B. E. Warm-water coral reefs and climate change. Science (80-. ). 350, 769–771 (2015).

9.
Rogers, A., Blanchard, J. L. & Mumby, P. J. Fisheries productivity under progressive coral reef degradation. J. Appl. Ecol. 55, 1041–1049 (2017).
Article  Google Scholar 

10.
Weijerman, M. et al. Evaluating management strategies to optimise coral reef ecosystem services. J. Appl. Ecol. 55, 1823–1833 (2017).
Article  Google Scholar 

11.
Jackson, J., Donovan, M., Cramer, K., Lam, V. & (editors). Status and Trends of Caribbean Coral Reefs: 1970–2012. (2014).

12.
Gardner, T. A., Cote, I. M., Gill, J. A., Grant, A. & Watkinson, A. R. Long-term region-wide declines in caribbean corals. Science (80-. ). 301, 958–961 (2003).

13.
Aronson, R., Bruckner, A., Moore, J., Precht, B. & Weil, E. Acropora palmata: The IUCN Red List of Threatened Species 2008: e.T133006A3536699. (2008). https://doi.org/10.2305/IUCN.UK.2008.RLTS.T133006A3536699.en.

14.
Aronson, R., Bruckner, A., Moore, J., Precht, B. & Weil, E. Acropora cervicornis. The IUCN Red List of Threatened Species 2008: e.T133381A3716457. (2008). https://doi.org/10.2305/IUCN.UK.2008.RLTS.T133381A3716457.en.

15.
Porto-Hannes, I. et al. Population structure of the corals Orbicella faveolata and Acropora palmata in the Mesoamerican Barrier Reef System with comparisons over Caribbean basin-wide spatial scale. Mar. Biol. 162, 81–98 (2015).
CAS  Article  Google Scholar 

16.
Keck, J., Houston, R. S., Purkis, S. & Riegl, B. M. Unexpectedly high cover of Acropora cervicornis on offshore reefs in Roatan (Honduras). Coral Reefs 24, 509 (2005).
ADS  Article  Google Scholar 

17.
Japaud, A., Bouchon, C., Manceau, J.-L. & Fauvelot, C. High clonality in Acropora palmata and Acropora cervicornis populations of Guadeloupe French Lesser Antilles. Mar. Freshw. Res. 66(847), 851 (2015).
Google Scholar 

18.
Baums, I. B., Miller, M. W. & Hellberg, M. E. Regionally isolated populations of an imperiled Caribbean coral Acropora palmata. Mol. Ecol. 14(1377), 1390 (2005).
Google Scholar 

19.
Lasker, H. R. & Coffroth, M. A. Responses of clonal reef taxa to environmental change. Am. Zool. 39, 92–103 (1999).
Article  Google Scholar 

20.
Honnay, O. & Bossuyt, B. Prolonged clonal growth: escape route or route to extinction?. Oikos 108, 427–432 (2005).
Article  Google Scholar 

21.
Teo, A. & Todd, P. A. Simulating the effects of colony density and intercolonial distance on fertilisation success in broadcast spawning scleractinian corals. Coral Reefs 37, 891–900 (2018).
ADS  Article  Google Scholar 

22.
Baums, I. B., Miller, M. W. & Hellberg, M. E. Geographic variation in clonal structure in a reef-building Caribbean coral Acropora palmata. Ecol. Monogr. 76(503), 519 (2006).
Google Scholar 

23.
Drury, C., Paris, C. B., Vassiliki, H. K. & Lirman, D. Dispersal capacity and genetic relatedness in Acropora cervicornis on the Florida Reef Tract. Coral Reefs 37, 585–596 (2018).
ADS  Article  Google Scholar 

24.
Reusch, T. B. H., Ehlers, A., Hämmerli, A. & Worm, B. Ecosystem recovery after climatic extremes enhanced by genotypic diversity. Proc. Natl. Acad. Sci. USA 102, 2826–2831 (2005).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

25.
Booth, R. E. & Grime, J. P. Effects of genetic impoverishment on plant community diversity. J. Ecol. 91, 721–730 (2003).
Article  Google Scholar 

26.
Drury, C., Greer, J. B., Baums, I., Gintert, B. & Lirman, D. Clonal diversity impacts coral cover in Acropora cervicornis thickets: Potential relationships between density, growth, and polymorphisms. Ecol. Evol. 9, 4518–4531 (2019).
PubMed  PubMed Central  Article  Google Scholar 

27.
Neigel, J. E. & Avise, J. C. Clonal diversity and population structure in a reef-building coral, acropora cervicornis: self-recognition analysis and demographic interpretation. Evolution (N. Y). 37, 437–453 (1983).

28.
Hughes, T. P. et al. Coral reefs in the anthropocene. Nature 546, 82–90 (2017).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

29.
Evensen, N., Doropoulos, C., Morrow, K., Motti, C. & Mumby, P. Inhibition of coral settlement at multiple spatial scales by a pervasive algal competitor. Mar. Ecol. Prog. Ser. 612, 29–42 (2019).
ADS  CAS  Article  Google Scholar 

30.
Goreau, T. J. & Hilbertz, W. Marine ecosystems restoration: costs and benefits for coral reefs. World Resour. Rev. 17, 375–409 (2005).
Google Scholar 

31.
Vollmer, S. V. & Kline, D. I. Natural disease resistance in threatened staghorn corals. PLoS ONE 3, 1–5 (2008).
Article  CAS  Google Scholar 

32.
Baums, I. B. A restoration genetics guide for coral reef conservation. Mol. Ecol. 17, 2796–2811 (2008).
PubMed  Article  PubMed Central  Google Scholar 

33.
Schopmeyer, S. A. et al. In situ coral nurseries serve as genetic repositories for coral reef restoration after an extreme cold-water event. Restor. Ecol. 20, 696–703 (2011).
Article  Google Scholar 

34.
Young, C. N., Schopmeyer, S. A. & Lirman, D. A review of reef restoration and coral propogation using the threatened genus Acropora in the Caribbean and Western Atlantic. Bull. Mar. Sci. 88, 1075–1098 (2012).
Article  Google Scholar 

35.
ICF. Instituto de Conservacion Forestal. Declaracion de Banco Cordelia Sitio de Importancia para la Vida Silvestre. Gaceta No. 32,816, 10 Mayo Del 2012. Acuerdo No. 021–2012. (2012).

36.
Riegl, B., Purkis, S. J., Keck, J. & Rowlands, G. P. Monitored and modeled coral population dynamics and the refuge concept. Mar. Pollut. Bull. 58, 24–38 (2009).
CAS  PubMed  Article  PubMed Central  Google Scholar 

37.
R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing (2020). https://www.R-project.org/.

38.
Soong, K. & Lang, J. C. Reproductive integration in reef corals. Biol. Bull. 183, 418–431 (1992).
CAS  PubMed  Article  PubMed Central  Google Scholar 

39.
Baums, I. B., Hughes, C. R. & Hellberg, M. E. Mendelian microsatellite loci for the Caribbean coral Acropora palmata. Mar. Ecol. Ser. 288, 115–127 (2005).
CAS  Article  Google Scholar 

40.
Baums, I. B., Devlin-Durante, K., Brown, L. & Pinzon, J. H. Nine novel, ploymorphic microsatellite markers for the study of threatened caribbean acroporid corals. Mol. Ecol. Resour. 9, 1155–1158 (2009).
CAS  PubMed  Article  PubMed Central  Google Scholar 

41.
Alberto, F. MsatAllele-1.0: An R package to visualize the binning of microsatellite alleles. J. Hered. 100, 394–397 (2009).

42.
Meirmans, P. G. & Van Tienderen, P. H. GENOTYPE and GENODIVE: Two programs for the analysis of genetic diversity of asexual organisms. Mol. Ecol. Notes 4, 792–794 (2004).
Article  Google Scholar 

43.
Peakall, R. & Smouse, P. E. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics 28, 2537–2539 (2012).

44.
Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 155, 945–959 (2000).
CAS  PubMed  PubMed Central  Google Scholar 

45.
Earl, D. A. & vonHoldt, B. M. STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361 (2012).

46.
Oksanen, J. et al. The Vegan Package. (2007).

47.
Vollmer, S. V. & Palumbi, S. R. Restricted gene flow in the Caribbean staghorn coral Acropora cervicornis: implications for the recovery of endangered reefs. J. Hered. 98, 40–50 (2007).
CAS  PubMed  Article  PubMed Central  Google Scholar 

48.
Van Woesik, R., Lacharmoise, F. & Köksal, S. Annual cycles of solar insolation predict spawning times of Caribbean corals. Ecol. Lett. 9, 390–398 (2006).
PubMed  Article  PubMed Central  Google Scholar 

49.
Fogarty, N. D., Vollmer, S. V. & Levitan, D. R. Weak prezygotic isolating mechanisms in threatened caribbean Acropora corals. PLoS One 7, (2012).

50.
Rodríguez-Martínez, R. E., Banaszak, A. T., McField, M. D., Beltrán-Torres, A. U. & Álvarez-Filip, L. Assessment of Acropora palmata in the mesoamerican reef system. PLoS ONE 9, 1–7 (2014).
Google Scholar 

51.
Aronson, J. & Alexander, S. Ecosystem restoration is now a global priority: time to roll up our sleeves. Restor. Ecol. 21, 293–296 (2013).
Article  Google Scholar 

52.
Perring, M. P. et al. Advances in restoration ecology: rising to the challenges of the coming decades. Ecosphere 6, 480–493 (2015).
Article  Google Scholar 

53.
Lirman, D. & Schopmeyer, S. Ecological solutions to reef degradation: optimizing coral reef restoration in the Caribbean and Western Atlantic. PeerJ 4, e2597 (2016).
PubMed  PubMed Central  Article  Google Scholar 

54.
Boström-Einarsson, L. et al. Coral restoration—a systematic review of current methods, successes, failures and future directions. PLoS ONE 15, 1–24 (2020).
Article  CAS  Google Scholar 

55.
Mijangos, J. L., Pacioni, C., Spencer, P. B. S. & Craig, M. D. Contribution of genetics to ecological restoration. Mol. Ecol. 24, 22–37 (2015).
PubMed  Article  PubMed Central  Google Scholar 

56.
Ladd, M. C., Miller, M. W., Hunt, J. H., Sharp, W. C. & Burkepile, D. E. Harnessing ecological processes to facilitate coral restoration. Front. Ecol. Environ. 16, 239–247 (2018).
Article  Google Scholar 

57.
Granado, R., Neta, L. C. P., Nunes-Freitas, A. F., Voloch, C. M. & Lira, C. F. Assessing genetic diversity after mangrove restoration in Brazil: Why is it so important? Diversity 10, (2018).

58.
Johnson, M. E. et al. Caribbean Acropora Restoration Guide: Best Practices for Propagation and Population Enhancement. (2011).

59.
Bland, L. M. et al. Using multiple lines of evidence to assess the risk of ecosystem collapse. Proc. R. Soc. B Biol. Sci. 284, 1–10 (2017).
Google Scholar 

60.
ICF. Plan de Gestión para el Manejo del Sitio de Importancia para la Vida Silvestre Banco Cordelia el Parque nacional Marino de Islas de la Bahia. Inst. Nac. Conserv. y Desarro. For. Areas Protegidas y Vida Silv. 111 (2013).

61.
Mcleod, E. et al. The future of resilience-based management in coral reef ecosystems. J. Environ. Manage. 233, 291–301 (2019).
PubMed  Article  PubMed Central  Google Scholar 

62.
Crouzeilles, R. et al. Ecological restoration success is higher for natural regeneration than for active restoration in tropical forests. Sci. Adv. 3, 1–8 (2017).
Article  Google Scholar  More

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

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

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

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

Full size image

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

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

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

Full size image

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

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

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

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

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

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

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

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

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

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

Full size image

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Full size image

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

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

Full size image

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