More stories

  • in

    Ongoing ecological and evolutionary consequences by the presence of transgenes in a wild cotton population

    In this study, we showed that the expression of cry and cp4-epsps genes in wild cotton altered the secretion of EFN, the associations with different ant species, and the levels of herbivore damage on target plants. Wcry constantly maintained a high production of EFN, regardless of the MeJA treatment, but nectar production was minimal in Wcp4-epsps. These changes in nectar inducibility seem to modify the composition of ant communities, foster the dominance of the generalist and defensive species C. planatus in Bt plants and the presence of ants without defensive role, M. ebeninum, in the herbicide tolerant genotype, while W plants had both defending species (C. planatus, C. rectangularis aulicus and P. gracilis) and invasive ant species (P. longicornis) in the same proportion. Furthermore, herbivore damage and its associated ant community were different according to the introgressed transgene.
    Wild and introgressed cotton do not display phenotypic equivalence in natural conditions
    In general, it has been assumed that introgressed and wild genotypes should display similar phenotypes in the absence of the selection agents targeted by transgenes. However, when we compared the control group and the three genotypes, we registered different nectar secretion patterns among them (Fig. 1). Similar results have been registered in populations of bt rice and glyphosate-tolerant sunflowers living in natural conditions where introgressed plants are different from their wild relatives5.
    Transgene expression modified indirect induced defences in wild cotton
    Most plants are able to induce responses after herbivore damage and/or phytohormone exogenous application (i.e. jasmonic acid, JA; methyl jasmonate, MeJA; and salicylic acid, SA)11,28,29. However, unlike wild plants without transgenes, individuals with transgenes were not sensitive to the induction treatment with MeJA for increasing their EFN production (Fig. 1). These results contrast with previous reports on cultivated varieties, such as Bt and glyphosate-resistant (cp4-epsps), in which direct defences such as gossypol terpenoids (160%), hemigossypolone (160%), helicoids 1|4 (213%) and indirect defenses, such as volatile compounds (VOCs) (171.2%) and extrafloral nectar (EFN) (133%), were reported to increase in plants sprinkled with JA and MeJA21,28,29,30.
    The inability of plants with transgenes to have the production of extrafloral nectar induced in them was related to different processes dependent on the identity of the transgenes in question. Whereas Wcry control plants had a high EFN production equivalent to the induced state of W plants, EFN production in Wcp4-epsps plants was inhibited. Contrasting these findings with results obtained under controlled conditions (i.e. greenhouse and crop conditions)3,21, we suggest that EFN production is linked to genotypes with transgenes and abiotic stress in the coastal dunes, because transgenes are connected to main metabolic pathways that respond to stressful conditions21.
    Wild cotton with cp4-epsps
    In the absence of herbicides acting as a selection agent, wild plants with cp4-epsps exhibited large differences compared to wild plants without them. Their low nectar production ( > 8 µg/mL) (Fig. 1) could be linked to the crosstalk between the jasmonate and the salicylate (SA) pathways (Fig. 4, orange and purple section). In G. hirsutum and other species, SA signalling has been proven to negatively affect JA signalling (e.g. Zea mays, Solanum lycopersicum, Nicotiana tabacum and Arabidopsis thaliana)31,32,33: therefore, we suggest an interference between the SA and JA pathways given previous reports that an over-expression of the cp4-epsps gene modifies the second part of the shikimate pathway (post-chorismate), which leads to the synthesis of essential amino acids as phenylalanine, tryptophan, or tyrosine, the latter being a precursor of benzoic acid BE, and SA34,35 (Fig. 4, purple section). This evidence highlights that hidden crosstalk effects among different metabolic pathways can scale up and modify plant phenotypes (e.g. extrafloral nectar production).
    Figure 4

    A diagram illustrating how the expression of cry (A) and cp4-epsps (B) in absence of their selection agent (pests and glyphosate) can affect the extrafloral nectar production. The extrafloral nectar (EFN) production is an induced defence that can be triggered by foliar herbivory, mechanical damage, and exogenous application of phytohormones (i.e. jasmonic acid, methyl jasmonate, and salicylic acid). These factors activate the octadecanoid pathway, and therefore, the production of extrafloral nectar, (A) aqua rectangle. The (C) section is an example of this reaction in a wild cotton plant (without transgenes). After damage, the key genes (yellow mesh) of the octadecanoid pathway are activated and produce extrafloral nectar. Another scenario is when the wild cotton expresses cry genes (A section), in this case, the key genes of the octadecanoid pathway interact synergistically with the cry transgene (green mesh). This triggers an over-expression of the production of EFN (aqua thick arrow), switching from inducible to constitutive responses. When the plants express cp4-epsps (B section), the production of extrafloral nectar is reduced or inhibited. A possible answer is an over-expression of the epsps gene (gold curve arrow), that increased production of salicylic acid which creates a crosstalk between shikimate and octadecanoid pathways (black cross-talk arrow). When the shikimate pathway is activated, the principal inducible defence is the production of volatile organic compounds (VOCs) (pink rectangle).

    Full size image

    Wild cotton with cry
    Wild cotton plants with cry genes continuously produced EFN as a constitutive defence (Fig. 1), in equivalent quantities as the induced state of W plants. EFN production is regulated by the octadecanoid signalling pathway, which can be activated by herbivore damage, mechanical damage, and phytohormones, such as JA and MeJA21,28 (Fig. 4, green section). However, for cotton, a specific elicitor is not necessary36. Four key genes for the synthesis of JA and MeJA have been described: AOS, AOC, HPL, and COI137. In Bt maize, studies comparing GM corn and its isogenic lines report an increase of 24% in phenols and 63% of DIMBOA (2,4-dihidroxi-7-metoxi-1,4-benzoxazin-3-ona; natural defences against lepidopteran herbivores)11. This is consistent with observations of a synergy between maize direct defences and Bt genes, after exogenous applications of JA (Fig. 4, orange section). Considering the latter, we suggest that Wcry cotton may present a similar response, as the genes activating the JA pathway are GhAOS and GhCOI1 (homologs to maize JA biosynthesis genes: ZmAOS and ZmCOI1), in addition to Ghppo1, which confers natural resistance to lepidopteran pest, such as H. armigera38. The interaction of cry with other genes could modify the production of EFN in Wcry plants.
    Effect of the transgenes’ expression on ants associated to wild cotton
    We identified eight species of ants harvesting EFN (Table 2), but with distinctive communities as a function of the plant genotype. This result suggests that the change in quantity, and possibly the composition and quality of EFN, can influence the ant community associated with G. hirsutum39,40,41.
    Changes in plant reward production could potentially compromise the attraction of natural enemies of herbivores42. In our study, the availability of EFN was modified. Although species richness was the same as in W plants (Table 2), the most abundant ant species associated with Wcp4-epsps plants, M. ebeninum, is considered a generalist species. Moreover, due to the lack of aggressive behaviour, this species does not represent an effective biotic defence43. The high abundance of this non-defensive species could be associated with the greater herbivore damage observed in Wcp4-epsps plants (Fig. 2). In contrast, W or Wcry plants showed a greater abundance of more aggressive ant species such as C. planatus, C. rectangulatus, and P. brunneus and significantly less herbivore damage.
    In Wcry cotton, the community of patrolling ants was mainly dominated by C. planatus, in both treatments (control and induction). Interestingly, although the amount of nectar did not vary between treatments, the abundance of ants was significantly different. The dominance of a single ant species could have benefited the plants with increased indirect defence, reducing herbivore damage and promoting a greater seed production per plant, as described in Turnera ulmifolia44, Schomburgkia tibicinis45, and Opuntia stricta42. However, considering the aggressive and dominant behaviour of C. planatus, there may be ecological costs through antagonistic relationships with pollinators. Ants can interrupt pollination and affect plant fitness25,46,47. The outcome of these mutualistic and antagonistic interactions requires further study.
    Effects of transgenes on herbivore damage
    Considering that the type of mutualism that cotton sustains with ants is defensive, we suggest that the change we observed in the composition of ants is likely to have influenced herbivore damage in the different genotypes, which in turn has the potential to reduce fitness as shown by other studies of cotton48,49,50. However, a study carried out on wild upland cotton reported that plants tolerate intermediate levels of leaf damage inflicted by leaf-chewing insects ( More

  • in

    Understanding the seasonality of performance resilience to climate volatility in Mediterranean dairy sheep

    1.
    Henry, B. K., Eckard, R. J. & Beauchemin, K. A. Review: Adaptation of ruminant livestock production systems to climate changes. Animal 12, s445–s456. https://doi.org/10.1017/S1751731118001301 (2018).
    CAS  Article  PubMed  Google Scholar 
    2.
    Sejian, V., Kumar, D., Gaughan, J. B. & Naqvi, S. M. K. Effect of multiple environmental stressors on the adaptive capability of Malpura rams based on physiological responses in a semi-arid tropical environment. J. Vet. Behav. 17, 6–13. https://doi.org/10.1016/j.jveb.2016.10.009 (2017).
    Article  Google Scholar 

    3.
    Segnalini, M., Bernabucci, U., Vitali, A., Nardone, A. & Lacetera, N. Temperature humidity index scenarios in the Mediterranean basin. Int. J. Biometeorol. 57, 451–458. https://doi.org/10.1007/s00484-012-0571-5 (2013).
    ADS  CAS  Article  PubMed  Google Scholar 

    4.
    Carabano, M. J., Bachagha, K., Ramon, M. & Diaz, C. Modeling heat stress effect on Holstein cows under hot and dry conditions: Selection tools. J. Dairy Sci. 97, 7889–7904. https://doi.org/10.3168/jds.2014-8023 (2014).
    CAS  Article  PubMed  Google Scholar 

    5.
    Finocchiaro, R., van Kaam, J. B., Portolano, B. & Misztal, I. Effect of heat stress on production of Mediterranean dairy sheep. J. Dairy Sci. 88, 1855–1864. https://doi.org/10.3168/jds.S0022-0302(05)72860-5 (2005).
    CAS  Article  PubMed  Google Scholar 

    6.
    Sevi, A. & Caroprese, M. Impact of heat stress on milk production, immunity and udder health in sheep: A critical review. Small Ruminant Res. 107, 1–7. https://doi.org/10.1016/j.smallrumres.2012.07.012 (2012).
    Article  Google Scholar 

    7.
    Giorgi, F. & Lionello, P. Climate change projections for the Mediterranean region. Glob. Planet. Change 63, 90–104. https://doi.org/10.1016/j.gloplacha.2007.09.005 (2008).
    ADS  Article  Google Scholar 

    8.
    Ozturk, T., Ceber, Z. P., Türkeş, M. & Kurnaz, M. L. Projections of climate change in the Mediterranean Basin by using downscaled global climate model outputs. Int. J. Climatol. 35, 4276–4292. https://doi.org/10.1002/joc.4285 (2015).
    Article  Google Scholar 

    9.
    Williams, C. M. et al. Understanding evolutionary impacts of seasonality: An introduction to the symposium. Integr. Comp. Biol. 57, 921–933. https://doi.org/10.1093/icb/icx122 (2017).
    Article  PubMed  PubMed Central  Google Scholar 

    10.
    Barash, H., Silanikove, N., Shamay, A. & Ezra, E. Interrelationships among ambient temperature, day length, and milk yield in dairy cows under a Mediterranean climate. J. Dairy Sci. 84, 2314–2320. https://doi.org/10.3168/jds.S0022-0302(01)74679-6 (2001).
    CAS  Article  PubMed  Google Scholar 

    11.
    Bernabucci, U. et al. Metabolic and hormonal acclimation to heat stress in domesticated ruminants. Animal 4, 1167–1183. https://doi.org/10.1017/S175173111000090X (2010).
    CAS  Article  PubMed  Google Scholar 

    12.
    Santana, M. L. Jr. et al. Detrimental effect of selection for milk yield on genetic tolerance to heat stress in purebred Zebu cattle: Genetic parameters and trends. J. Dairy Sci. 98, 9035–9043. https://doi.org/10.3168/jds.2015-9817 (2015).
    CAS  Article  PubMed  Google Scholar 

    13.
    Escarcha, J., Lassa, J. & Zander, K. Livestock under climate change: A systematic review of impacts and adaptation. Climate https://doi.org/10.3390/cli6030054 (2018).
    Article  Google Scholar 

    14.
    van der Werf, J., Graser, H.-U., Frankham, R. & Gondro, C. Adaptation and Fitness in Animal Populations: Evolutionary and Breeding Perspectives on Genetic Resource Management (Springer, New York, 2009).
    Google Scholar 

    15.
    van der Waaij, E. H. A resource allocation model describing consequences of artificial selection under metabolic stress. J. Anim. Sci. 82, 973–981. https://doi.org/10.2527/2004.824973x (2004).
    Article  PubMed  Google Scholar 

    16.
    Colditz, I. G. & Hine, B. C. Resilience in farm animals: Biology, management, breeding and implications for animal welfare. Anim. Prod. Sci. https://doi.org/10.1071/an15297 (2016).
    Article  Google Scholar 

    17.
    Berghof, T. V. L., Poppe, M. & Mulder, H. A. Opportunities to improve resilience in animal breeding programs. Front. Genet. 9, 692. https://doi.org/10.3389/fgene.2018.00692 (2019).
    Article  PubMed  PubMed Central  Google Scholar 

    18.
    Mulder, H. A. Genomic selection improves response to selection in resilience by exploiting genotype by environment interactions. Front. Genet. 7, 178. https://doi.org/10.3389/fgene.2016.00178 (2016).
    CAS  Article  PubMed  PubMed Central  Google Scholar 

    19.
    Bishop, S. C. A consideration of resistance and tolerance for ruminant nematode infections. Front. Genet. 3, 168. https://doi.org/10.3389/fgene.2012.00168 (2012).
    Article  PubMed  PubMed Central  Google Scholar 

    20.
    Knap, P. Breeding robust pigs. Aust. J. Exp. Agric. https://doi.org/10.1071/EA05041 (2005).
    Article  Google Scholar 

    21.
    Marjanovic, J., Mulder, H. A., Ronnegard, L. & Bijma, P. Modelling the co-evolution of indirect genetic effects and inherited variability. Heredity 121, 631–647. https://doi.org/10.1038/s41437-018-0068-z (2018).
    Article  PubMed  PubMed Central  Google Scholar 

    22.
    Mulder, H. A., Bijma, P. & Hill, W. G. Prediction of breeding values and selection responses with genetic heterogeneity of environmental variance. Genetics 175, 1895–1910. https://doi.org/10.1534/genetics.106.063743 (2007).
    CAS  Article  PubMed  PubMed Central  Google Scholar 

    23.
    Mulder, H. A., Ronnegard, L., Fikse, W. F., Veerkamp, R. F. & Strandberg, E. Estimation of genetic variance for macro- and micro-environmental sensitivity using double hierarchical generalized linear models. Genet. Sel. Evol. 45, 23. https://doi.org/10.1186/1297-9686-45-23 (2013).
    Article  PubMed  PubMed Central  Google Scholar 

    24.
    Westneat, D. F., Wright, J. & Dingemanse, N. J. The biology hidden inside residual within-individual phenotypic variation. Biol. Rev. Camb. Philos. Soc. 90, 729–743. https://doi.org/10.1111/brv.12131 (2015).
    Article  PubMed  Google Scholar 

    25.
    Strandberg, E., Brotherstone, S., Wall, E. & Coffey, M. P. Genotype by environment interaction for first-lactation female fertility traits in UK dairy cattle. J. Dairy Sci. 92, 3437–3446. https://doi.org/10.3168/jds.2008-1844 (2009).
    CAS  Article  PubMed  Google Scholar 

    26.
    Bernabucci, U. et al. The effects of heat stress in Italian Holstein dairy cattle. J. Dairy Sci. 97, 471–486. https://doi.org/10.3168/jds.2013-6611 (2014).
    CAS  Article  PubMed  Google Scholar 

    27.
    Carabaño, M. J. et al. Breeding and genetics symposium: Breeding for resilience to heat stress effects in dairy ruminants. A comprehensive review. J. Anim. Sci. https://doi.org/10.2527/jas2016.1114 (2017).
    Article  PubMed  Google Scholar 

    28.
    Ramon, M., Diaz, C., Perez-Guzman, M. D. & Carabano, M. J. Effect of exposure to adverse climatic conditions on production in Manchega dairy sheep. J. Dairy Sci. 99, 5764–5779. https://doi.org/10.3168/jds.2016-10909 (2016).
    CAS  Article  PubMed  Google Scholar 

    29.
    Sanchez-Molano, E. et al. Genetic analysis of novel phenotypes for farm animal resilience to weather variability. BMC Genet. 20, 84. https://doi.org/10.1186/s12863-019-0787-z (2019).
    Article  PubMed  PubMed Central  Google Scholar 

    30.
    Knap, P. W. & Su, G. Genotype by environment interaction for litter size in pigs as quantified by reaction norms analysis. Animal 2, 1742–1747. https://doi.org/10.1017/S1751731108003145 (2008).
    CAS  Article  PubMed  Google Scholar 

    31.
    Kolmodin, R., Strandberg, E., Madsen, P., Jensen, J. & Jorjani, H. Genotype by environment interaction in nordic dairy cattle studied using reaction norms. Acta Agric. Scand. Sect. A Anim. Sci. 52, 11–24. https://doi.org/10.1080/09064700252806380 (2002).
    Article  Google Scholar 

    32.
    Rauw, W. M. & Gomez-Raya, L. Genotype by environment interaction and breeding for robustness in livestock. Front. Genet. 6, 310. https://doi.org/10.3389/fgene.2015.00310 (2015).
    CAS  Article  PubMed  PubMed Central  Google Scholar 

    33.
    Stinchcombe, J. R., Function-valued Traits Working Group, G. & Kirkpatrick, M. Genetics and evolution of function-valued traits: Understanding environmentally responsive phenotypes. Trends Ecol. Evol. 27, 637–647. https://doi.org/10.1016/j.tree.2012.07.002 (2012).
    Article  PubMed  Google Scholar 

    34.
    Reeves, J. L. et al. Seasonal temperature and precipitation effects on cow–calf production in northern mixed-grass prairie. Livestock Sci. 155, 355–363. https://doi.org/10.1016/j.livsci.2013.04.015 (2013).
    Article  Google Scholar 

    35.
    Lateef, A., Das, H., Panchasara, H. H., Nilufar, H. & Sanap, M. J. Seasonal effects on milk yield, erythrocytic and leukocytic indices of Kankrej cattle (Bos indicus). Vet. World 7, 472–477. https://doi.org/10.14202/vetworld.2014.472-477 (2014).
    Article  Google Scholar 

    36.
    Angilletta, M. J. Thermal adaptation: A theoretical and empirical synthesis. Therm. Adapt. Theor. Empir. Synth. https://doi.org/10.1093/acprof:oso/9780198570875.001.1 (2009).
    Article  Google Scholar 

    37.
    Kingsolver, J., Diamond, S. & Gomulkiewicz, R. Integrative Organismal Biology (eds. Martin, L. B., Ghalambor, C. K. & Woods, H. A.) 39–53 (Wiley, New York, 2015).

    38.
    Hoffmann, I. Climate change and the characterization, breeding and conservation of animal genetic resources. Anim. Genet. 41(Suppl 1), 32–46. https://doi.org/10.1111/j.1365-2052.2010.02043.x (2010).
    Article  PubMed  Google Scholar 

    39.
    Molik, E. et al. Effect of day length and exogenous melatonin on chemical composition of sheep milk. Arch. Anim. Breed. 54, 177–187. https://doi.org/10.5194/aab-54-177-2011 (2011).
    CAS  Article  Google Scholar 

    40.
    Molik, E., Misztal, T. & Dorota A. The Effect of Physiological and Environmental Factors on the Prolactin Profile in Seasonally Breeding Animals. https://www.intechopen.com/books/prolactin/the-effect-of-physiological-and-environmental-factors-on-the-prolactin-profile-in-seasonally-breedin (2013).

    41.
    Peana, I., Fois, G. & Cannas, A. Effects of heat stress and diet on milk production and feed and energy intake of Sarda ewes. Italian J. Anim. Sci. https://doi.org/10.4081/ijas.2007.1s.577 (2010).
    Article  Google Scholar 

    42.
    Peana, I. et al. Cold markedly influences milk yield of Sardinian dairy sheep farms. Italian J. Anim. Sci. 6, 580–580. https://doi.org/10.4081/ijas.2007.1s.580 (2007).
    Article  Google Scholar 

    43.
    do Amaral, B. C. et al. Heat-stress abatement during the dry period: Does cooling improve transition into lactation?. J. Dairy Sci. 92, 5988–5999. https://doi.org/10.3168/jds.2009-2343 (2009).
    CAS  Article  PubMed  Google Scholar 

    44.
    do Amaral, B. C. et al. Heat stress abatement during the dry period influences prolactin signaling in lymphocytes. Domest. Anim. Endocrinol. 38, 38–45. https://doi.org/10.1016/j.domaniend.2009.07.005 (2010).
    CAS  Article  PubMed  Google Scholar 

    45.
    do Amaral, B. C. et al. Heat stress abatement during the dry period influences metabolic gene expression and improves immune status in the transition period of dairy cows. J. Dairy Sci. 94, 86–96. https://doi.org/10.3168/jds.2009-3004 (2011).
    CAS  Article  PubMed  Google Scholar 

    46.
    Serradilla, J. M. et al. Characterisation of Goats’ Response to Heat Stress: Tools to Improve Heat Tolerance. https://www.intechopen.com/books/goat-science/characterisation-of-goats-response-to-heat-stress-tools-to-improve-heat-tolerance (2018).

    47.
    Beilharz, R. G., Luxford, B. G. & Wilkinson, J. L. Quantitative genetics and evolution: Is our understanding of genetics sufficient to explain evolution?. J. Anim. Breed. Genet. 110, 161–170. https://doi.org/10.1111/j.1439-0388.1993.tb00728.x (1993).
    CAS  Article  PubMed  Google Scholar 

    48.
    Rauw, W. M. “Introduction” in Resource Allocation Theory Applied to Farm Animal Production (ed. Rauw, W.M.) 1–21 (CAB International Publishing, Wallingford, 2008).

    49.
    Garner, J. B. et al. Genomic selection improves heat tolerance in dairy cattle. Sci. Rep. 6, 34114. https://doi.org/10.1038/srep34114 (2016).
    ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

    50.
    International Committee for Animal Recording (ICAR). ICAR Guidelines: Dairy Sheep and Goats. http://www.icar.org (2017).

    51.
    Basdagianni, Z., Sinapis, E. & Banos, G. Evaluation of reference lactation length in Chios dairy sheep. Animal 13(1), 1–7. https://doi.org/10.1017/S1751731118000769 (2018).
    Article  PubMed  Google Scholar 

    52.
    Misztal, I. et al. BLUPF90 and related programs (BGF90). In Proceedings of 7th World Congress on Genetics Applied to Livestock Production, Vol. 743 (2002).

    53.
    R Core Team. R: A language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, 2019). http://www.r-project.org/index.html.
    Google Scholar 

    54.
    International Committee for Animal Recording (ICAR). ICAR Guidelines: Dairy Sheep and Goats. http://www.icar.org (2020)

    55.
    Gilmour, A. R., Gogel, B. J., Cullis, B. R., Welham, S. J. & Thompson, R. ASReml User Guide Release 4.1 Structural Specification. (2015). More

  • in

    The rise of angiosperms strengthened fire feedbacks and improved the regulation of atmospheric oxygen

    1.
    Bond, W. J., Woodward, F. I. & Midgley, G. F. The global distribution of ecosystems in a world without fire. N. Phytol. 165, 525–538 (2005).
    CAS  Article  Google Scholar 
    2.
    Krawchuk, M. A., Moritz, M. A., Parisien, M.-A., Van Dorn, J. & Hayhoe, K. Global pyrogeography: the current and future distribution of wildfire. PLoS ONE https://doi.org/10.1371/journal.pone.0005102 (2009).

    3.
    Pausas, J. G. & Keeley, J. E. A burning story: the role of fire in the history of life. Bio. Sci. 59, 593–601 (2009).
    Google Scholar 

    4.
    Pausas, J. G. & Ribeiro, E. The global fire-productivity relationship. Glob. Ecol. Biogeog. 22, 728–736 (2013).
    Article  Google Scholar 

    5.
    Archibald, S. et al. Biological and geophysical feedbacks with fire in the Earth system. Environ. Res. Lett. 13, 033003 (2018).
    ADS  Article  Google Scholar 

    6.
    Pausas, J. G., Bradstock, R. A., Keith, D. A. & Keeley, J. E. Plant functional traits in relation to fire in crown-fire ecosystems. Ecology 85, 1085–1100 (2004).
    Article  Google Scholar 

    7.
    Lenton, T. M. in Fire Phenomena and the Earth System: An Interdisciplinary Guide to Fire Science (ed Belcher, C. M.) 289–309 (John Wiley and Sons, Chichester, UK, 2013).

    8.
    van der Werf, G. R. et al. Global fire emissions and the contribution of deforestation, savanna, forest, agricultural, and peat fires 1997−2009. Atmos. Chem. Phys. Discuss. 10, 11707–11735 (2010).
    ADS  Article  CAS  Google Scholar 

    9.
    Santín, C. et al. Towards a global assessment of pyrogenic carbon from vegetation fires Glob. Chang. Biol. 22, 76–91 (2016).
    Article  Google Scholar 

    10.
    Bond, W. J. & Keeley, J. E. Fire as a global herbivore: the ecology and evolution of flammable ecosystems. Trends Ecol. Evol. 20, 387–394 (2005).
    PubMed  Article  Google Scholar 

    11.
    Belcher, C. M., Collinson, M. E. & Scott, A. C. in C. M. Belcher (ed) Fire Phenomena and the Earth System: an Interdisciplinary Guide to Fire Science 229–249 (Wiley, Oxford, 2016).

    12.
    Keeley, J. E., Pausas, J. G., Rundel, P. W., Bond, W. J. & Bradstock, R. A. Fire as an evolutionary pressure shaping plant traits. Trends Plant Sci. 16, 406–411 (2011).
    CAS  PubMed  Article  Google Scholar 

    13.
    He, T., Pausas, J. G., Belcher, C. M., Schwilk, D. W. & Lamont, B. B. Fire-adapted traits of Pinus arose in the fiery Cretaceous. N. Phytol. 194, 751–759 (2012).
    Article  Google Scholar 

    14.
    Belcher, C. M. & Hudspith, V. A. Changes to cretaceous surface fire behaviour influenced the spread of the early angiosperms. N. Phytol. 213, 1521–1532 (2016).
    Article  CAS  Google Scholar 

    15.
    Bond, W. J. & Scott, A. C. Fire and the spread of flowering plants in the Cretaceous. N. Phytol. 188, 1137–1150 (2010).
    Article  Google Scholar 

    16.
    Keeley, J. E. & Rundel, P. W. Fire and the Miocene expansion of C4 grasslands. Ecol. Lett. 8, 683–690 (2005).
    Article  Google Scholar 

    17.
    Belcher, C. M. & McElwain, J. C. Limits for combustion in low O2 redefine paleoatmospheric predictions for the Mesozoic. Science 321, 1197–1200 (2008).
    ADS  CAS  PubMed  Article  Google Scholar 

    18.
    Belcher, C. M., Yearsley, J. M., Hadden, R. M., McElwain, J. C. & Rein, G. Baseline intrinsic flammability of Earth’s ecosystems estimated from paleoatmospheric oxygen over the past 350 million years. Proc. Natl Acad. Sci. USA 107, 22448–22453 (2010).
    ADS  CAS  PubMed  Article  Google Scholar 

    19.
    Holland, H. D. The oxygenation of the atmosphere and oceans. Philos. Trans. Roy. Soc. B 361, 903–915 (2006).
    CAS  Article  Google Scholar 

    20.
    Goldblatt, C., Lenton, T. M. & Watson, A. J. Bistability of atmospheric oxygen and the Great Oxidation. Nature 443, 683–686 (2006).
    ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

    21.
    Alcott, L. J., Mills, B. J. W. & Poulton, S. W. Stepwise Earth oxygenation is an inherent property of global biogeochemical cycling. Science 366, 1333–1337 (2019).
    ADS  CAS  PubMed  Article  Google Scholar 

    22.
    Lenton, T. M. & Watson, A. J. Biotic enhancement of weathering, atmospheric oxygen and carbon dioxide in the Neoproterozoic. Geophys. Res. Lett. 31, L05202 (2004).
    ADS  Article  CAS  Google Scholar 

    23.
    Lenton, T. M. et al. Earliest land plants created modern levels of atmospheric oxygen. Proc. Natl Acad. Sci. USA 113, 9704–9709 (2016).
    ADS  CAS  PubMed  Article  Google Scholar 

    24.
    Lenton, T. M. & Watson, A. J. Redfield revisited 2. What regulates the oxygen content of the atmosphere? Glob. Biogeochem. Cycles, B 14, 249–268 (2000).
    ADS  CAS  Article  Google Scholar 

    25.
    Watson, A. J. Consequences for the biosphere of forest and grassland fires. PhD thesis (University of Reading, Reading, UK, 1978).

    26.
    Kump, L. R. Terrestrial feedback in atmospheric oxygen regulation by fire and phosphorus. Nature 335, 152–154 (1988).
    ADS  CAS  Article  Google Scholar 

    27.
    Burdige, D. J. Burial of terrestrial organic matter in marine sediments: a reassessment. Glob. Biogeochem. Cyc. 19, (2005).

    28.
    Van Cappellen, P. & Ingall, E. D. Redox stabilisation of the atmosphere and oceans by phosphorus-limited marine productivity. Science 271, 493–496 (1996).
    ADS  PubMed  Article  Google Scholar 

    29.
    Kump, L. R. & Mackenzie, F. T. Regulation of atmospheric O2: feedback in the microbial feedbag. Science 271, 459–460 (1996).
    ADS  CAS  Article  Google Scholar 

    30.
    Canfield, D. E. Sulfur isotopes in coal constrain the evolution of the Phaneroiz sulphur cycle. Proc. Natl Acad. Sci. USA 110, 8443–8446 (2013).
    ADS  CAS  PubMed  Article  Google Scholar 

    31.
    Walker, J. C. G. Stability of atmospheric oxygen. Am. J. Sci. 274, 193–214 (1974).
    ADS  CAS  Article  Google Scholar 

    32.
    Lasaga, A. C. & Ohmoto, H. The oxygen geochemical cycle: dynamics and stability. Geochim. Cosmochim. Acta 66, 361–381 (2002).
    ADS  CAS  Article  Google Scholar 

    33.
    Holland, H. D. The Chemical Evolution of the Atmosphere and Oceans: Princeton Series in Geochemistry (Princeton University Press, Princeton, 1984).

    34.
    Berner, R. A. Atmospheric oxygen over Phanerozoic time. Proc. Natl Acad. Sci. USA 96, 10955–10957 (1999).
    ADS  CAS  PubMed  Article  Google Scholar 

    35.
    Lovelock, J. Gaia: a new look at life on Earth (Oxford University Press, New York, 1979).

    36.
    Mahowald, N. M. et al. Impact of biomass burning emissions and land use change on Amazonian atmospheric cycling and deposition of phosphorus. Glob. Biogeochem. Cycles 19, 1–15 (2005).
    Google Scholar 

    37.
    Bergman, N. M., Lenton, T. M. & Watson, A. J. COPSE: a new model of biogeochemical cycling over Phanerozoic time. Am. J. Sci. 304, 397–437 (2004).
    ADS  CAS  Article  Google Scholar 

    38.
    Robinson, J. M. Phanerozoic O2 variation, fire, and terrestrial ecology. Palaeogeogr. Palaeoclimatol. Palaeoecol. 75, 223–240 (1989).
    Article  Google Scholar 

    39.
    Midgley, J. J. & Bond, W. J. in Fire Phenomena and the Earth System: An Interdisciplinary Guide to Fire Science (ed Belcher C. M.) 125–134 (Chichester, UK: John Wiley and Sons, 2013).

    40.
    Lenton, T. M., Daines, S. J. & Mills, B. J. W. COPSE reloaded: an improved model of biogeochemical cycling over Phanerozoic time. Earth Sci. Rev. 178, 1–28 (2018).
    ADS  CAS  Article  Google Scholar 

    41.
    Mills, B. J. W. et al. Modeliing the long-term carbon cycle, atmospheric CO2, and Earth surface temperature from late Neoproterozoic to present day. Gond. Res. 67, 172–186 (2019).
    ADS  CAS  Article  Google Scholar 

    42.
    Berner, R. A. Biogeochemical cycles of carbon and sulphur and their effect on atmospheric oxygen over Phanerozoic time. Glob. Planet. Chang. 75, 97–122 (1989).
    ADS  Article  Google Scholar 

    43.
    Berner, R. A. Modeling atmospheric O2 over Phaneroic time. Geochim. Cosmochim. Acta 65, 685–694 (2001).
    ADS  CAS  Article  Google Scholar 

    44.
    Lenton, T. M. The role of land plants, phosphorus weathering and fire in the rise and regulation of atmospheric oxygen. Glob., Chang. Biol. 7, 613–629 (2001).
    ADS  Article  Google Scholar 

    45.
    Berner, R. A. Phanerozoic atmospheric oxygen: new results using the GEOCARBSULF model. Am. J. Sci. 289, 333–361 (2009).
    ADS  Article  Google Scholar 

    46.
    Berner, R. A. & Canfield, D. E. A new model for atmospheric oxygen over Phanerozoic time. Am. J. Sci. 289, 333–361 (1989).
    ADS  CAS  PubMed  Article  Google Scholar 

    47.
    Glasspool, I. J. & Scott, A. C. Phanerozoic concentrations of atmospheric oxygen reconstructed from sedimentary charcoal. Nat. Geosci. 3, 627–630 (2010).
    ADS  CAS  Article  Google Scholar 

    48.
    Mills, B. J. W., Belcher, C. M., Lenton, T. M. & Newton, R. J. A modeling case for high atmospheric oxygen concentrations during the Mesozoic and Cenozoic. Geology 44, 1023–1026 (2016).
    ADS  CAS  Article  Google Scholar 

    49.
    Crane, P. R. & lidgard, S. Angiosperm diversification and paleolatitudinal gradients in Cretaceous floristic diversity. Science 246, 675–678 (1989).
    ADS  CAS  PubMed  Article  Google Scholar 

    50.
    Lidgard, S. & Crane, P. R. Angiosperm diversification and Cretaceous floristic trends – a comparison of palynofloras and leaf floras. Paleobiology 16, 77–93 (1990).
    Article  Google Scholar 

    51.
    Friis, E. M., Crane, P. R. & Pederson, J. R. Early Flowers and Angiosperm Evolution (Cambridge University Press, New York, 2011).

    52.
    Bond, W. J. The tortoise and the hare: ecology of angiosperm dominance and gymnosperm persistence. Biol. J. Linn. Soc. 36, 227–249 (1989).
    Article  Google Scholar 

    53.
    Brodribb, T. J., Feild, T. S. & Jordan, G. J. Leaf maximum photosynthetic rate and venation are linked by hydraulics. Plant Physiol. 144, 1890–1898 (2007).
    CAS  PubMed  PubMed Central  Article  Google Scholar 

    54.
    Field, T. S. et al. Fossil evidence for Cretaceous escalation un angiosperm leaf vein evolution. Proc. Natl Acad. Sci. USA 108, 863–8366 (2011).
    Google Scholar 

    55.
    Field, T. S., Arens, N. C., Doyle, J. A. & Dawson, T. E. Dark and disturbed: a new image of early angiosperm ecology. Paleobiology 30, 82–107 (2004).
    Article  Google Scholar 

    56.
    Mills, B. et al. Changing tectonic controls on the long-term carbon cycle from Mesozoic to present. Geochem. Geophys. Geosyst. 15, 4866–4884 (2014).
    ADS  CAS  Article  Google Scholar 

    57.
    Schneider, H. et al. Ferns diversified in the shadow of angiosperms. Nature 428, 553–557 (2004).
    ADS  CAS  PubMed  Article  Google Scholar 

    58.
    Brodribb, T. J. & Field, T. S. Leaf hydraulic evolution led a surge in leaf photosynthetic capacity during early angiosperm diversification. Ecol. Lett. 13, 175–183 (2010).
    PubMed  Article  Google Scholar 

    59.
    Royer, D. L., Miller, I. M., Peppe, D. J. & Hickey, L. J. Leaf economic traits from fossils support a weedy habit for early angiosperms. Am. J. Bot. 97, 438–445 (2010).
    PubMed  Article  Google Scholar 

    60.
    Stevens P. F. Angiosperm phylogeny website. Version 12, [Online]. http://www.mobot.org/MOBOT/research/APweb/ (2012).

    61.
    Schwilk, D. W. & Ackerly, D. D. Flammability and serotiny as strategies: correlated evolution in pines. Oikos 94, 326–336 (2001).
    Article  Google Scholar 

    62.
    Lamont, B. B. & Downes, K. S. Fire-stimulated flowering among resprouters and geophytes in Australia and South Africa. Pl Ecol. 212, 2111–2125 (2011).
    Article  Google Scholar 

    63.
    Janssen, T. & Bremer, K. The age of major monocot groups inferred from 800+rbcL sequences. Bot. J. Linn. Soc. 146, 385–398 (2004).
    Article  Google Scholar 

    64.
    Krassilov, V. & Voynets, Y. Weedy Albian angiosperms. Acta Palaeobot. 48, 151–169 (2008).
    Google Scholar 

    65.
    Lamont, B. B. & He, T. Fire adapted Gondanan angiosperm floras evolved in the Cretaceous. Bmc. Evol. Biol. 12, 223 (2012).
    PubMed  PubMed Central  Article  Google Scholar 

    66.
    Bowman, D. M. J. S., French, B. J. & Prior, L. D. Have plants evolved to self-immolate? Front. Plant. Sci. https://doi.org/10.3389/fpls.2014.00590 (2014).

    67.
    Crisp, M. D., Burrows, G. E., Cook, l. G., Thornhill, A. H. & Bowman, D. M. J. S. Flammable biomes dominated by eucalypts originated at the Cretaceous-Palaeogene boundary. Nat. Commun 2, 1–8 (2011).
    Article  CAS  Google Scholar 

    68.
    He, T., Lamont, B. B. & Downes, K. S. Banksia born to burn. N. Phyt. 191, 184–196 (2011).
    Article  Google Scholar 

    69.
    Watson, A. J. & Lovelock, J. E. in Fire Phenomena and the Earth System: An Interdisciplinary Guide to Fire Science (ed Belcher C. M.) 273–287 (John Wiley and Sons, Chichester, UK, 2013).

    70.
    Tiffney, B. H. Seed size, dispersal syndromes, and the rise of the angiosperms: evidence and hypothesis. Ann. Miss. Bot. Gard. 71, 551–576 (1984).
    Article  Google Scholar 

    71.
    Wheeler, E. A. & Baas, P. A survey of the fossil record for dicotlydonous wood and its significant for evolutionary and ecological wood anatomy. Int. Assoc. Wood. Anatom. Bull. 12, 271–332 (1991).
    Google Scholar 

    72.
    Wing, S. L. & Boucher, L. Ecological aspects of the Cretaceous flowering plant radiation. Ann. Rev. Earth Plan. Sci. 26, 379–421 (1998).
    ADS  CAS  Article  Google Scholar 

    73.
    Johnson, K. R. & Ellis, B. A tropical rainforest in Colorado, 1.4 million years after the Cretaceous-Tertiary boundary. Science 296, 2379–2383 (2002).
    ADS  CAS  PubMed  Article  Google Scholar 

    74.
    Spicer, R. A., McA, Rees, P. & Chapman, J. L. Cretaceous phytogeography and climate signals. Roy. Soc. Proc. B 341, 277–286 (1993).
    Google Scholar 

    75.
    Beerling, D. J. & Woodward, F. I. Vegetation and the Terrestrial Carbon Cycle: Modelling the First 400 Million Years. (Cambridge, Cambridge and New York, 2001).

    76.
    Bond, W. J. & Midgley, J. J. Fire and the angiosperm revolutions. Int. J. Plant Sci. 173, 1–16 (2012).
    Article  Google Scholar 

    77.
    Wing, S. L. et al. Late Paleocene fossils from the Cerrejón formation, Colombia, are the earliest record of Neotropical rainforest. Proc. Natl Acad. Sci. USA 106, 18627–18632.

    78.
    Boyce, C. K., Brodribb, T., Feild, T. S. & Zweiniecki, M. J. Angiosperm leaf vein evolution was physiologically and environmentally transformative. Proc. Roy. Soc. B 276, 1771–1776 (2009).
    Article  Google Scholar 

    79.
    Balch, J. K. and five others. Human-started wildfires expand the fire niche across the United States. Proc. Natl Acad. Sci. USA 114, 2946–2951 (2017).
    ADS  CAS  PubMed  Article  Google Scholar 

    80.
    Donovan, G. H. & Brown, T. C. Be careful what you wish for: The legacy of Smokey Bear. Front. Ecol. Environ. 5, 73–79 (2007).
    Article  Google Scholar  More

  • in

    Susceptibility of Pimephales promelas and Carassius auratus to a strain of koi herpesvirus isolated from wild Cyprinus carpio in North America

    Collection of wild carp from a CyHV-3-exposed population
    This study was carried out in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All protocols for sampling, procedures and experimental endpoints involving live fish conducted in this study were approved by the Institutional Animal Care & Use Committee (IACUC), University of Minnesota (St. Paul, Minnesota, USA), under the approval numbers IACUC-1806-36036A and 1808-36276A. Experiments were performed in compliance with the ARRIVE guidelines on animal research32.
    Wild carp were sampled from Lake Elysian (Waseca County, Minnesota, Coordinates: 44.178144, − 93.69066) by boat electrofishing from September 3 to 9, 2019 (Fig. 1a). This lake was expected to have a CyHV-3-exposed carp population following a confirmed outbreak in 20173. Captured wild adult carp (n = 116) were euthanized by immersion in a solution of ~ 3 mL/L pure clove oil (90% Eugenol; Velona, Elk Grove Village, IL, USA) for 15 min and transported on ice to the University of Minnesota for necropsy. Brain, gill and kidney tissues from up to three carp were pooled in a 1:5 (weight:volume) dilution of Hank’s Balanced Salt Solution (HBSS; Cellgro, Lincoln, NE, USA) containing 100 IU/mL of Penicillin and Streptomycin and maintained at a pH of 7.4 at 4 °C for 24 h prior to preparation for qPCR and cell culture screening for CyHV-3 (described below). Gill tissues from ten freshly-dead carp obtained from a shallow bay in the Southern portion of the lake were also obtained and pooled by five individuals for a total of two sample pools.
    Figure 1

    (a) Generated using ArcMap (v10.8.1, https://desktop.arcgis.com/en/arcmap/), shows the approximate locations of sampling effort and mortality observations on Lake Elysian. Bathymetric contours indicate depth in 5 ft increments. (b, c) Pathology of a representative individual wild carp sampled from Lake Elysian. Arrows on (b, c) denote frayed fins (vermillion), loss of mucosal layer (white), loss of scales and epidermis (black), enopthalmia (bluish green), gill necrosis (sky blue).

    Full size image

    An additional 17 wild carp collected as part of the previously described sampling event were placed in an aerated live well and transferred to the Minnesota Aquatic Invasive Species Research Center’s Containment Laboratory (MCL). These carp were housed in a ~ 1400 L tank with flow through well water at 20 °C and treated with 0.6% aquarium salt once per day. Carp were acclimated for 1 day and then anesthetized via immersion in a solution of 100 µL/L of clove oil and uniquely marked using colored injectable elastomer (Northwest Marine Technology, Anacortes, WA, USA). Additionally, a small portion (~ 0.2 cm2) of each carp’s gills were sampled for qPCR screening for CyHV-3 and tested immediately. Carp determined to be CyHV-3 negative (n = 12) were euthanized following testing. Carp determined to be CyHV-3-positive (n = 5) by specific qPCR were held for a total of 5 days, during which, water temperature was gradually increased to 26 °C in order to increase viral shedding. CyHV-3-positive carp gill biopsies were again sampled and screened on the fifth day to identify carp with high qPCR copy numbers. All CyHV-3-positive carp were then euthanized, and the brain, gill and kidney tissues were removed as previously described. Pooled tissues from two wild carp with clinical signs consistent with KHVD (Fig. 1b,c) and with high qPCR copy numbers, were subjected to cell culture immediately following necropsy. In addition, a 10 g portion of this pooled tissue was processed and used to challenge naive carp in the in-vivo infection model. Tissues were homogenized in a 1:5 volume of HBSS containing 100 IU/mL Penicillin and Streptomycin (pH = 7.4). The sample was centrifuged at 2360 × g at 25 °C for 10 min, then the supernatant was passed through a 0.45 µm syringe filter.
    In-vivo infection trial
    To increase the potential of obtaining an isolate of CyHV-3, naïve carp previously determined to be CyHV-3 negative by qPCR, were challenged with CyHV-3-positive tissue homogenates obtained from wild carp. Two naïve carp, purchased from Osage Catfisheries (Osage Beach, MO, USA), were pair housed in a 60 L aquarium with flow through well water (flow rate = 3–4 tank volumes/h) at 21–22 °C. Aquaria were set up with a standpipe drain covered by a cylindrical wire screen filter of approximately 15 cm in length and 4.4 cm in diameter. Additionally, a PVC pipe section of 15 cm in length and 10 cm in diameter was added to each tank for shelter. Each carp was exposed to 0.5 mL of CyHV-3-positive tissue homogenate by IP-injection and monitored for signs of disease for 6 days and then euthanized. Pooled samples of brain, gill and kidney tissue were subjected to qPCR and cell culture analysis. Following cell culture analysis (below) a second infection trial was performed to satisfy River’s postulates (i.e. that CyHV-3 isolated from wild diseased carp would cause similar disease in naïve carp)33. Two additional naïve carp purchased from Osage Catfisheries were IP-injected with 0.5 mL of CyHV-3-positive (qPCR and cell culture positive) cell culture supernatant. Carp were housed and observed for disease signs as previously described for 11 days and then sacrificed. Pooled samples of brain, gill, and kidney then were tested by CyHV-3-specific qPCR to confirm the presence of CyHV-3.
    Cell culture analysis
    CCB cells were maintained in Eagle’s Minimum Essential Medium (EMEM) containing Eagles’s salts (Sigma, St. Louis, MO, USA), 10% fetal bovine serum (FBS), 1% non-essential amino acids (NEAA, Sigma), 2 mM l-glutamine and glucose (Sigma) up to 4.5 g/L. The KF-1 cells were cultured in EMEM containing Eagles’s salts (Sigma), 10% FBS and 0.025 M HEPES. Penicillin 100 U/L and streptomycin 0.1 mg/L (Sigma) were used as an anti-bacterial agent in both cell culture media and the cells were maintained at 25 °C.
    Cell culture methods to isolate CyHV-3 were performed according to the US Fish and Wildlife Service and American Fisheries Society-Fish Health Section Blue Book34. Briefly, pooled tissues were homogenized in Hank’s Balanced Salt Solution (HBSS; Cellgro) and centrifuged at 2360 × g for 15 min. The inoculum was added to the 24-well plates with 80% confluent cell cultures in two dilutions, (1/10 and 1/100) and incubated at 25 °C for 14 days. A blind passage was performed for an additional 14 days if no cytopathic effects (CPE) were observed on the first passage. If CPE was observed during the first passage, then the second passage was performed in a 25 cm2 flask. The virus was harvested when CPE reached 70–80% of the monolayer. The infected cultures were exposed to two freeze/thaw cycles at − 80 °C, and then centrifuged at 3800 × g for 15 min at 4 °C. The clarified supernatants and pellets were collected and stored at − 80 °C.
    Whole-genome sequencing and sequence analysis
    Whole-genome sequencing was performed at the University of Minnesota Veterinary Diagnostic Laboratory for genetic characterization of the CyHV-3 isolate (KHV/Elysian/2019) obtained from wild carp. In brief, after CCB cells, infected with wild carp tissues, reached 80% CPE, the supernatant was collected and stored at − 80 °C. The frozen supernatant was freeze-thawed three times, and centrifuged at 2896 × g for 25 min at 4 °C. Nucleic acid purification of CCB cell culture supernatant was done using a QIAamp MinElute Virus Spin Kit (Qiagen, Hilden, Germany) following manufacturer instructions. The extracted nucleic acids were subjected to library preparation using Nextera Flex DNA library kit (Illumina, San Diego, CA, USA) following manufacturer instructions. The library was normalized according to the median fragment size measured by Tape Station 2.0 (Agilent, Santa Clara, CA, USA) and library concentration measured by Qubit. The library was submitted to the University of Minnesota Genomic Center (UMGC) for sequencing via MiSeq V3 (2X75-bp) paired end chemistry.
    Raw fastq files were trimmed to remove Illumina adapters using Trimmomatic (v 0.39) with a minimum quality score of 20. Then, bowtie2 (v 2.3.5) was used to remove host contamination and unmapped reads were used for assembly with SPAdes (v3.13.0) with k-mer values of 29, 33 and 55 with the options “careful with a minimum coverage of 5 reads per contig”. Then contigs were searched into the RefSeq viral and non-redundant protein reference databases using Diamond BLASTx with an e-value of 1e − 5 for significant hits. Taxon assignments were made with MEGAN6 Community Edition according to the lowest-common-ancestor algorithm. ORFs prediction and genome annotation were done using Prokka (v1.14.5). The resulting alignment (GenBank accession no. MT914509) was aligned with 19 other CyHV-3 genomes available on NCBI using Mafft (v7) with the FFT-NS-2 alignment strategy and the following parameters: scoring matrix BLOUSUM62, gap open penalty 1.53, offset value 0. Model selection, maximum likelihood (ML) tree construction, and calculation of bootstrap values were done with R 4.0 (R Software) using phangorn (v2.5.5). ML trees were constructed using the top scoring model (GTR + G + I) and 100 bootstrap replicates were generated to assess the reliability of clades obtained in the tree. Additionally, this genome assembly was compared with the previously reported thymidine kinase gene sequence obtained from carp sampled during a large mortality event in Lake Elysian in 2017 (F36, GenBank accession no. MK987089).
    Investigation of species specificity
    To investigate the host range of KHV/Elysian/2019, six carp purchased from Osage Catfisheries, previously determined to be CyHV-3-negative by qPCR, were intraperitoneally (IP) injected with 0.5 mL of the filtered tissue homogenate material (Fig. 2a). The IP-injected carp (IP-carp) were housed as previously described for 9 days prior to their use in the cohabitation trial (Fig. 2b). The IP-carp were monitored twice daily for signs of disease. After 9 days the gills, skin and vent of each IP-carp was swabbed aseptically with a single sterile cotton swab (Dynarex, Orangeburg, NY, USA) for determination of viral load by qPCR. FHM and goldfish were challenged with CyHV-3 via cohabitation. One cohabitation tank (tank A) contained ten naïve FHM, five naïve sentinel carp (S-carp) and three IP-carp (Fig. 2a). One cohabitation tank (tank B) contained ten naïve goldfish, five naive S-carp and three of the IP-carp. S-carp were included in each tank setup to act as a positive control for within-tank transmission of CyHV-3. Two additional negative control tanks with the same stocking density and conditions contained ten naïve FHM (tank D) and ten naïve goldfish (tank E), as well as eight naïve carp (confirmed to be CyHV-3-negative by specific qPCR). Average standard length and weight for fishes used in these experiments was 13 cm and 64 g for carp, 7 cm and 13 g for FHM, and 10 cm and 38 g for goldfish. All tanks consisted of ~ 60 L aquaria with flow-through well water as previously described. Fishes were fed a commercial feed (Skretting classic trout, Skretting, Tooele, UT, USA) daily and monitored twice daily to observe changes to fish health. IP-carp that died during the trial were allowed to remain in the tank for 24 h prior to removal for necropsy, but any morbidity or mortality of other experimental groups were immediately removed and necropsied.
    Figure 2

    (a) Shows a schematic of the cohabitation disease trial. Vermillion arrows denote inoculation of IP-carp with CyHV-3 positive tissue homogenate, blue arrows denote introduction of IP carp for cohabitation with fishes in experimental tanks, and the reddish purple arrow indicates the tissue origin of CyHV-3-positive S-carp. (b) Shows a schematic of experimental flow through chambers with black arrows indicating the direction of water flow. (c) Shows a time-line of various samples.

    Full size image

    At 0, 3, 6, 9, 12, and 15 days post exposure (dpe) by cohabitation, five FHM, five goldfish, and all IP-carp and S-carp from each tank were anesthetized by immersion in a buffered solution of 100 mg/L of MS-222 and the gills, skin and vent of each fish was swabbed with a sterile swab for determination of viral load by qPCR (Fig. 2c). For FHM and goldfish, the five individuals were randomly sampled at each time-point. Additionally, the wire screen filter of the outflow standpipe was swabbed at the same intervals during the course of the trial to evaluate loading of CyHV-3 DNA in the environment. All swabs were stored at − 20 °C in individual plastic bags until nucleic acid extraction could be performed. At 11 dpe, half of the FHM and goldfish from cohabitation tanks were euthanized by immersion in a buffered solution of 3 g/L of MS-222 and necropsied (Fig. 2c). The remaining FHM and goldfish were maintained until 20 dpe and then euthanized and necropsied. To visually record the presence of gross pathology, representative IP carp, and fish from cohabitation groups (S-carp, FHM, and goldfish) were randomly selected and photographed at 0 and 6 dpe in a small glass aquarium (Fig. 3).
    Figure 3

    Representative fishes photographed before and after exposure to CyHV-3. Note, fishes photographed at 0 dpe may not be the same individual as those at 6 dpe. dpe days post exposure via cohabitation, IP-carp intraperitoneally injected carp, S-carp cohabitated sentinel carp, FHM fathead minnow. Arrows denote frayed fins (vermillion), loss of mucosal layer (white), scale pocket edema (black). Additionally, normal morphological features of mature male fathead minnows are indicated for nuptial tubercles (bluish green), and nape pads varying in prominence (reddish purple).

    Full size image

    For each necropsied fish, wet mounts of gill and skin scrapes were viewed at 40× magnification to identify potential parasitic infections. Then the skin of each necropsied fish was rinsed briefly with 70% ETOH and clean water. Brain, gill, kidney and skin tissue were collected individually for each fish and split into two duplicate samples. The first sample duplicates were placed in Whirl–Pak sample bags (Nasco, Fort Atkinson, WI, USA) and preserved at − 20 °C until nucleic acid extraction and screening for CyHV-3 DNA was performed. The second sample duplicates were placed in 1 mL of RNAlater solution (Ambion) in 1.5 mL microcentrifuge tubes (Globe Scientific, Mahwah, NJ, USA) and frozen at − 20 °C. An individual FHM and goldfish from each time-point (11- and 20-dpe) was preserved in 10% NBF (TissuePro, Gainesville, FL, USA) for histological analysis. Individual representatives of each species from control tanks and moribund S-carps from each experimental tank were also preserved for histological analysis.
    Due to the detection of CyHV-3 DNA in a single FHM in tank A, a second trial with FHM (tank C) was performed as described previously (Fig. 2a). Brain, gill, kidney, and skin tissue from two S-carp exposed in the first trial with disease signs and positive qPCR test for CyHV-3 (tank A) were pooled, homogenized and filtered as previously described. Three new carp purchased from Osage Catfisheries were IP injected with 0.5 mL of this tissue homogenate and maintained as previously described for 9 days prior to screening for CyHV-3 by qPCR and used in the cohabitation trial. All other conditions and procedures were done as described for the first cohabitation trial with the following exceptions. In the second trial, portions of brain, gill, kidney and skin tissues obtained from a moribund S-carp at 5 dpe and four FHM at 11 dpe, respectively, were pooled as previously described and subjected to cell culture. Additionally, duplicate swabs from the tank C outflow standpipe filter were obtained and preserved in 1 mL of RNAlater solution (Sigma) as previously described for tissue samples.
    Nucleic acid purification using chelex resin and detection of CyHV-3 by qPCR
    For nucleic acid purification, chelex resin (Sigma) was used as described by Zida et al.35 and briefly summarized here. For pooled tissue samples, approximately 100 mg of each tissue was homogenized in 1 mL of nuclease free water (NFW) and then centrifuged, with 50 μL of the resulting supernatant later used as starting material. For swabs, the cotton end was cut off and vortexed, then centrifuged and finally the cotton was removed leaving the supernatant. For each sample type, 150 μL of chilled 80% ETOH was added, then centrifuged and the supernatant removed. Samples were allowed to air dry for 10 min to remove residual ETOH. 150 μL of 20% Chelex was added to each sample and vortexed. Samples were then incubated at 90 °C for 20 min and centrifuged and immediately used for qPCR.
    A Taqman probe-based qPCR was used for the detection of CyHV-3 DNA targeting the ORF89 gene36 using a StepOnePlus thermocycler with default settings (Applied Biosystems). Nucleic acid purifications from all samples were screened for CyHV-3 using a PrimeTime gene expression master mix kit (Integrated DNA Technologies, Coralville, IA, USA), with each reaction containing 400 nM of primers (KHV-86f: GAC-GCC-GGA-GAC-CTT-GTG, KHV-163r: CGG-GTT-GTT-ATT-TTT-GTC-CTT-GTT) and 250 nM of the probe (KHV-109p: [TAMRA] CTT-CCT-CTG-CTC-GGC-GAG-CAC-G-[IBRQ]. The reaction mix was subjected to an initial denaturation at 95 °C for 3 min, followed by 40 cycles of denaturation at 95 °C for five sec and annealing at 60 °C for 30 s. A threshold cycle of 38 was used as a cut off. The standard curve for quantification of CyHV-3 genomes was performed using a laboratory synthesized DNA fragment containing the ORF89 sequence as previously described by Padhi et al.3. The results for virus load are presented as the number of viral copies per mL of tissue supernatant. All samples obtained from FHM and goldfish were tested in triplicate with the exception of samples that had positive qPCR Ct values, which were re-tested up to six times.
    RNA purification and reverse transcription polymerase chain reaction (RT-PCR)
    Individual tissues of preserved brain, gill, kidney, and skin from one representative S-carp from each experimental tank (A, B and C) were selected as positive controls for CyHV-3 mRNA detection (total of 12 tissue samples). All preserved tissue samples from FHM or goldfish which had at least one positive qPCR test were also screened for CyHV-3 mRNA to determine if replicating virus was present (total of eight tissue samples). Additionally, preserved swabs of the outflow standpipe filter were also screened. For RNA purification, RNA was extracted from tissues using the RNeasy Mini Kit (Qiagen) according to the manufacturer instructions for animal tissues, using ~ 30 mg tissue samples preserved in RNAlater. For swabs, cotton was cut from the end of the swab and used as the starting material. CyHV-3 mRNA was detected using the RT-PCR developed by Yuasa et al.29 with the primers, (KHV RT F3: GCC-ATC-GAC-ATC-ATG-GTG-CA, KHV RT R1: AAT-GCC-GCT-GGA-AGC-CAG-GT). The RT-PCR was performed using a One-step RT-PCR kit (Qiagen) according to the manufacturer instructions. The reaction mix was subjected to a single step of reverse transcription at 50 °C for 30 min and denaturation at 95 °C for 15 min, followed by 40 cycles of: 94 °C for 30 s, 65 °C for 30 s, 72 °C for one minute and a final extension step was 72 °C for 10 min. PCR products were separated and visualized on 2% agarose gels containing 0.75 μg/mL ethidium bromide (Genesee Scientific, San Diego, CA, USA). PCR products for carp, FHM and goldfish templates (clear band at the 219 bp location) were cut from gels and purified by precipitation with a 20% PEG, 2.5 M NaCl solution. Purified RT-PCR products were subjected to Sanger sequencing at the University of Minnesota Genomics Center (UMGC). Sequences were trimmed and analyzed using 4 peaks (v1.8) and consensus sequences were generated using BioEdit (v7.2.1). Sequence identities were compared with available reference sequences by BLASTn analysis of the National Center of Biotechnology sequence database.
    Histology
    Histology was used to demonstrate the presence or absence of lesions in cohabitation disease trial specimens. Histological samples of gill tissue were prepared from formalin-fixed samples of representative fishes of each species from trial and control tanks. Gill samples were dissected from formalin-fixed specimens and decalcified in 10% ethylenediaminetetraacetic acid (EDTA) for 10 days. Following decalcification, samples were dehydrated in an ethanol series to 100% ethanol, infiltrated with toluene, and subsequently embedded in paraffin. Paraffin sections were cut at 6 µm thickness using a Leica Jung 820 Histocut Rotary Microtome and mounted on slides. Sections were stained with Hematoxylin and Eosin using a protocol modified from Humasson37.
    Statistical analysis
    R 4.0 (R Software) was used for statistical analysis and data presentation. CyHV-3 qPCR copy numbers are presented as averages of all positive tests for samples with duplicate tests and were Log transformed prior to statistical testing. Significant differences (p  More

  • in

    Mating and starvation modulate feeding and host-seeking responses in female bed bugs, Cimex lectularius

    The bed bug Cimex lectularius is an obligate ectoparasite that engages in recurrent blood-feeding events throughout its lifetime, punctuated by sheltering some distance from the host. This complex lifestyle requires coordination of diverse on-host and off-host behavioral repertoires, including host-seeking, blood-feeding, mating, oviposition, and aggregation to sustain development, reproduction, and survival11,13. Adult female bed bugs are expected to monitor their changing physiological state and nutritional and reproductive needs, as well as environmental cues such as host availability, and express specific behaviors in an adaptive coordinated manner that supports reproductive success while minimizing fitness costs. While developing and validating a vertical Y-olfactometer for bed bugs, we observed a mating-dependent behavioral shift where 47% of mated females responded to human skin odor, but none of 20 unmated females responded to the same olfactory stimuli23. In the present report, we followed-up on this observation because, to our knowledge, no experimental studies have investigated whether bed bugs modulate host-seeking and blood-feeding behaviors with changes in their mating and nutritional states.
    Consistent with our previous results 23, we found that 64–69% of mated females responded to human skin odor 8 days after their last blood meal (Fig. 3). In contrast, only 24% of unmated females responded, revealing again that host-seeking behavior is profoundly modulated by the mating state of the bed bug. We observed higher host-seeking response rates in this study than in the previous report23, possibly because we used slightly different durations of starvation and olfactometer conditions.
    We speculated that the behavioral differences between unmated and mated females were related to the rate of processing of the blood meal, which would be affected in turn by the rate of oocyte maturation and oviposition. Specifically, we hypothesized that because unmated females resorb their eggs24, they have lower nutritional requirements and reduced metabolic rate25, and therefore engage in less host-seeking and blood-feeding. To test this hypothesis, we extended the starvation period for both mated and unmated females. We found that the length of starvation had different effects on the host-seeking and feeding responses of mated and unmated females. Whereas host-seeking gradually increased in unmated females (24 to 60%) through 40 days of starvation, mated females became less responsive to host cues with longer starvation. These observations were consistent with the blood-feeding assays—with longer starvation, more unmated females fed and they took larger blood meals, even though they did not fully process some of their previous blood meal, as evidenced by their unfed body mass (Fig. 1b).
    These results are consistent with the hypothesis that the differences between mated and unmated females are driven by the accelerated egg maturation and oviposition cycle of mated females. Mated females initiated egg laying 4 days after ingesting a blood meal, oviposited on average 15.3 eggs per female, and ceased oviposition 6 days later, 10 days after ingesting a blood meal (Fig. 4). Thus, 8 days after a blood meal, mated females digested most of the blood (Fig. 1b), oviposited most of their eggs, and became highly motivated to host-seek to support a second oviposition cycle. Indeed, both laboratory and field observations showed that, given the opportunity, mated females accept a second blood meal while the first blood meal is still being digested and females feed every 2.5 days on average11,13,26. This is in contrast with other hematophagous insects, such as mosquitoes, where female host-seeking and feeding are suppressed for several days after a blood meal, until she completes laying a batch of eggs18,27.
    The strategy of unmated females was to feed little when the host is available, but with longer starvation periods, they became more stimulated to host-seek and ingest increasingly larger blood meals. This strategy is likely driven by much-reduced nutritional demands related to resorption of oocytes, which allow unmated females to digest the blood meal more slowly and use it for somatic maintenance rather than egg maturation. As stated by Davis in24: “If the female has fed but has not been inseminated, egg maturation will proceed to an early stage and then the yolk material will be resorbed; thus virgin females avoid the waste of producing unfertilized eggs.” Metabolic differences between mated and unmated females also support the idea that unmated females avoid wasting resources—unmated females have reduced metabolic rates after feeding compared with mated females25,28. Overall, this strategy would result in fewer host-seeking forays and less blood ingested by unmated females (Fig. 1b,c), until they mate and are stimulated to mature and oviposit fertile eggs. A similar strategy appears to operate in the closely related kissing bug R. prolixus, where virgin females remained unresponsive or even repelled by host-associated cues (CO2 and heat) until 20 days after ingesting a blood meal19. Rhodnius is a much larger hemipteran than C. lectularius, it takes larger blood meals, and likely requires more time to digest the blood.
    Surprisingly, host-seeking declined in mated females that were starved for 30 or 40 days, and this was especially apparent in females housed with fertile males (mated-long treatment) (Fig. 3). Two factors might account for this observation. The first is female aging and senescence, as 40 days of starvation in these females was beyond the 35-day median survival of females in this treatment group, and 100% of these females died by day 49 (Fig. 5). Thus, the females that we bioassayed 40 days after they ingested a blood meal were likely weak and less responsive to olfactory stimuli. This reduction in the host-seeking and blood-feeding responses of older mated females might be associated with their higher metabolic rate, senescence, or aging that could negatively affect olfactory responses, as shown in the D. melanogaster29.
    The second factor that likely underlies their early senescence is the unusual extra-genitalic, hemocoelic (traumatic) insemination in C. lectularius. Females housed with fertile males would receive multiple copulations that represent constant harassment, stress, and physical damage. These interactions with sexually aggressive males reduced their median adult lifespan by 63% (mated-long treatment) relative to females that were housed with fertile males for only 24 h and then with another female (mated-short treatment) (Fig. 5). These findings are consistent with previous studies in bed bugs on the adverse effects of multiple copulations on female lifespan13,14,17,30,31. We also observed that mated-long females assumed a “refusal” posture, protecting the ectospermalege from the males (Fig. 6a, Supplementary Video S1). It is possible, however, that males might circumvent these defensive postures by puncturing the intersegmental membranes away from the specialized ectospermalege, as shown in the closely related Cimex hemipterus32, and thus cause substantial damage to the female. This injurious effect of males on females might shape female behavioral responses in the field, but these responses were obviously constrained under our experimental conditions. For example, mated females might leave aggregations to avoid males, as demonstrated experimentally in laboratory assays33. Mated females might also seek refugia that are too narrow to accommodate themselves as well as courting males. Whether these evasive behaviors occur under field conditions will require further observations.
    We found no significant difference in female fecundity in the first oviposition cycle (~ 10 days post blood meal) between two treatment groups—females with long- and short-term presence of fertile males that represented high and low mating rates, respectively (Fig. 4). It is important to emphasize, however, that this experiment was limited to a single feeding and only one oviposition cycle. The lifespan of mated females was dramatically reduced by both the high and low mating rates, but significantly more so by high mating rate (Fig. 5), consistent with previous observations30. Morrow and Arnqvist30 also found that while elevated mating rate shortened female lifespan, it did not affect lifetime egg production. Although we used a different experimental design and did not measure lifetime fecundity, these findings suggest that mated females preferentially direct resources to egg maturation at the risk of significantly reduced lifespan, a strategy that requires close monitoring of physiological state and environmental conditions34.
    Remarkably, we also detected a significant effect on females of non-copulatory harassment by males. Females housed with a male that could not copulate (intromittent organ ablated) for only 24 h and then with another female (unmated-short treatment) lived to a median age of 111 days (100% dead by day 142), whereas females housed continuously with an infertile male (unmated-long treatment) lived to a median age of only 73 days (100% died by day 104). This 34% decline in expected lifespan, independent of copulation and egg production, can be attributed to male-specific harassment (Fig. 5). Males engage in a stereotyped behavior where the male repeatedly mounts the female’s dorsum, bends his abdomen to her ventral surface, and probes the female’s sternites with the paramere. We observed that 50% of the unmated-long and 62.5% of mated-long females exhibited a ‘refusal’ posture at least once in their lifetime thereby making the ectospermalege inaccessible to males, while none of the unmated-short or mated-short females displayed this behavior (Fig. 6b). This behavior may be similar to a behavior noted by N. Davis [in11, p. 171], but not described, that “starved females exhibit a distinct avoidance of mating”. The expression of this refusal behavior in virgin females that obviously need to mate is particularly surprising, suggesting that male harassment may be especially detrimental to the metabolic budget of starved females that feed less and endeavor to conserve energy. Unfortunately, our experimental design did not enable us to determine whether co-habitation with a female also would decrease survivorship of starved females relative to solitary females. It is possible that general disturbance of the starved female causes her to move and expend energy, which in turn reduces her lifespan. In this context, it is worth noting that by adding a female to the mated-short treatment, after the male was removed, the presence of the extra female might have decreased survivorship and thus resulted in underestimating the difference between the mated-long and mated-short treatments. More

  • in

    Family graveyards form underappreciated local plant diversity hotspots in China's agricultural landscapes

    1.
    Tilman, D., Isbell, F. & Cowles, J. M. Biodiversity and ecosystem functioning. Annu. Rev. Ecol. Evol. Syst. 45, 471–493. https://doi.org/10.1146/annurev-ecolsys-120213-091917 (2014).
    Article  Google Scholar 
    2.
    Isbell, F. et al. Biodiversity increases the resistance of ecosystem productivity to climate extremes. Nature 526, 574–577. https://doi.org/10.1038/nature15374 (2015).
    ADS  CAS  Article  PubMed  Google Scholar 

    3.
    Wood, S. A. et al. Functional traits in agriculture: agrobiodiversity and ecosystem services. Trends Ecol. Evol. 30, 531–539. https://doi.org/10.1016/j.tree.2015.06.013 (2015).
    Article  PubMed  Google Scholar 

    4.
    Salek, M. et al. Bringing diversity back to agriculture: smaller fields and non-crop elements enhance biodiversity in intensively managed arable farmlands. Ecol. Ind. 90, 65–73. https://doi.org/10.1016/j.ecolind.2018.03.001 (2018).
    Article  Google Scholar 

    5.
    Bianchi, F. J., Booij, C. J. & Tscharntke, T. Sustainable pest regulation in agricultural landscapes: a review on landscape composition, biodiversity and natural pest control. Proc. Biol. Sci. 273, 1715–1727. https://doi.org/10.1098/rspb.2006.3530 (2006).
    CAS  Article  PubMed  PubMed Central  Google Scholar 

    6.
    Newbold, T. et al. Has land use pushed terrestrial biodiversity beyond the planetary boundary? A global assessment. Science 353, 288–291 (2016).
    ADS  CAS  Article  Google Scholar 

    7.
    Foley, J. A. et al. Global consequences of land use. Science 309, 570–574. https://doi.org/10.1126/science.1111772 (2005).
    ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

    8.
    Green, R. E., Cornell, S. J., Scharlemann, J. P. W. & Balmford, A. Farming and the fate of wild nature. Science 307, 550–555. https://doi.org/10.1126/science.1106049 (2005).
    ADS  CAS  Article  PubMed  Google Scholar 

    9.
    Carvalheiro, L. G., Seymour, C. L., Nicolson, S. W., Veldtman, R. & Clough, Y. Creating patches of native flowers facilitates crop pollination in large agricultural fields: mango as a case study. J. Appl. Ecol. 49, 1373–1383. https://doi.org/10.1111/j.1365-2664.2012.02217.x (2012).
    Article  Google Scholar 

    10.
    Lindborg, R., Plue, J., Andersson, K. & Cousins, S. A. O. Function of small habitat elements for enhancing plant diversity in different agricultural landscapes. Biol. Conserv. 169, 206–213. https://doi.org/10.1016/j.biocon.2013.11.015 (2014).
    Article  Google Scholar 

    11.
    Knappova, J., Hemrova, L. & Muenzbergova, Z. Colonization of central European abandoned fields by dry grassland species depends on the species richness of the source habitats: a new approach for measuring habitat isolation. Landsc. Ecol. 27, 97–108. https://doi.org/10.1007/s10980-011-9680-5 (2012).
    Article  Google Scholar 

    12.
    Mendenhall, C. D., Karp, D. S., Meyer, C. F. J., Hadly, E. A. & Daily, G. C. Predicting biodiversity change and averting collapse in agricultural landscapes. Nature 509, 213. https://doi.org/10.1038/nature13139 (2014).
    ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

    13.
    Hunter, M. L. et al. Conserving small natural features with large ecological roles: a synthetic overview. Biol. Conserv. 211, 88–95. https://doi.org/10.1016/j.biocon.2016.12.020 (2017).
    Article  Google Scholar 

    14.
    Batáry, P. et al. Effect of conservation management on bees and insect-pollinated grassland plant communities in three European countries. Agric. Ecosyst. Environ. 136, 35–39. https://doi.org/10.1016/j.agee.2009.11.004 (2010).
    Article  Google Scholar 

    15.
    McGlaughlin, M. E. et al. Do the island biogeography predictions of MacArthur and Wilson hold when examining genetic diversity on the near mainland California Channel Islands? Examples from endemic Acmispon (Fabaceae). Bot. J. Linn. Soc. 174, 289–304. https://doi.org/10.1111/boj.12122 (2014).
    Article  Google Scholar 

    16.
    Wilcove, D. S., McLellan, C. H. & Dobson, A. P. Habitat Fragmentation in the Temperate Zone (Sinauer Associates Inc, Sunderland, Massachusetts, 1986).
    Google Scholar 

    17.
    Baur, B. & Erhardt, A. Habitat fragmentation and habitat alterations: principal threats to most animal and plant species. GAIA Ecol. Perspect. Sci. Soc. 4, 221–226 (1995).
    Google Scholar 

    18.
    Mac Arthur, R. H. & Wilson, E. O. The Theory of Island Biogeography. Monographs in Population Biology. (Princeton University Press, Princeton, N. J., 1967).

    19.
    Lindgren, J. P. & Cousins, S. A. O. Island biogeography theory outweighs habitat amount hypothesis in predicting plant species richness in small grassland remnants. Landsc. Ecol. 32, 1895–1906. https://doi.org/10.1007/s10980-017-0544-5 (2017).
    Article  Google Scholar 

    20.
    Li, L. Distribution pattern of plant diversity and vegetation construction in field margins and homegardens Doctor thesis, China Agriculture University, (2014).

    21.
    Li, P. et al. Possibilities and requirements for introducing agri-environment measures in land consolidation projects in China, evidence from ecosystem services and farmers’ attitudes. Sci. Total Environ. 650, 3145–3155. https://doi.org/10.1016/j.scitotenv.2018.10.051 (2019).
    ADS  CAS  Article  PubMed  Google Scholar 

    22.
    Haddad, N. M. et al. Plant species loss decreases arthropod diversity and shifts trophic structure. Ecol. Lett. 12, 1029–1039. https://doi.org/10.1111/j.1461-0248.2009.01356.x (2009).
    Article  PubMed  Google Scholar 

    23.
    Haddad, N. M., Tilman, D., Haarstad, J., Ritchie, M. & Knops, J. M. H. Contrasting effects of plant diversity and composition on insect communities: a field experiment. Am. Nat. 158, 17–35 (2001).
    CAS  Article  Google Scholar 

    24.
    Hertzog, L. R., Meyer, S. T., Weisser, W. W. & Ebeling, A. Experimental manipulation of grassland plant diversity induces complex shifts in aboveground arthropod diversity. PLoS ONE 11, e0148768. https://doi.org/10.1371/journal.pone.0148768 (2016).
    CAS  Article  PubMed  PubMed Central  Google Scholar 

    25.
    Southwood, T. R. E., Brown, V. K. & Reader, P. M. The relationships of plant and insect diversities in succession. Biol. J. Lin. Soc. 12, 327–348. https://doi.org/10.1111/j.1095-8312.1979.tb00063.x (1979).
    Article  Google Scholar 

    26.
    Hector, A. et al. Plant diversity and productivity experiments in European grasslands. Science 286, 1123–1127. https://doi.org/10.1126/science.286.5442.1123 (1999).
    CAS  Article  PubMed  Google Scholar 

    27.
    Tilman, D., Hill, J. & Lehman, C. Carbon-negative biofuels from low-input high-diversity grassland biomass. Science 314, 1598–1600. https://doi.org/10.1126/science.1133306 (2006).
    ADS  CAS  Article  PubMed  Google Scholar 

    28.
    Cardinale, B. J. et al. Impacts of plant diversity on biomass production increase through time because of species complementarity. Proc. Natl. Acad. Sci. U. S. A. 104, 18123–18128. https://doi.org/10.1073/pnas.0709069104 (2007).
    ADS  Article  PubMed  PubMed Central  Google Scholar 

    29.
    Hautier, Y. et al. Anthropogenic environmental changes affect ecosystem stability via biodiversity. Science 348, 336–340. https://doi.org/10.1126/science.aaa1788 (2015).
    ADS  CAS  Article  PubMed  Google Scholar 

    30.
    Tilman, D., Reich, P. B. & Knops, J. M. H. Biodiversity and ecosystem stability in a decade-long grassland experiment. Nature 441, 629–632. https://doi.org/10.1038/nature04742 (2006).
    ADS  CAS  Article  PubMed  Google Scholar 

    31.
    Klein, A. M. et al. Importance of pollinators in changing landscapes for world crops. Proc. Biol. Sci. 274, 303–313. https://doi.org/10.1098/rspb.2006.3721 (2007).
    Article  PubMed  Google Scholar 

    32.
    Hoehn, P., Tscharntke, T., Tylianakis, J. M. & Steffan-Dewenter, I. Functional group diversity of bee pollinators increases crop yield. Proc. R. Soc. B Biol. Sci. 275, 2283–2291. https://doi.org/10.1098/rspb.2008.0405 (2008).
    Article  Google Scholar 

    33.
    Kleijn, D. et al. Ecological effectiveness of agri-environment schemes in different agricultural landscapes in the Netherlands. Conserv. Biol. 18, 775–786. https://doi.org/10.1111/j.1523-1739.2004.00550.x (2004).
    Article  Google Scholar 

    34.
    Steffan-Dewenter, I. & Tscharntke, T. Succession of bee communities on fallows. Ecography 24, 83–93 (2001).
    Article  Google Scholar 

    35.
    Steffan-Dewenter, I., Klein, A.-M., Gaebele, V., Alfert, T. & Tscharntke, T. Bee Diversity and Plant-Pollinator Interactions in Fragmented Landscapes (UNIV Chicago Press, 2006).

    36.
    Bruun, H. H. Patterns of species richness in dry grassland patches in an agricultural landscape. Ecography 23, 641–650. https://doi.org/10.1034/j.1600-0587.2000.230601.x (2000).
    Article  Google Scholar 

    37.
    Öster, M., Cousins, S. A. O. & Eriksson, O. Size and heterogeneity rather than landscape context determine plant species richness in semi-natural grasslands. J. Veg. Sci. 18, 859–868 (2007).
    Article  Google Scholar 

    38.
    Tscharntke, T., Klein, A. M., Kruess, A., Steffan-Dewenter, I. & Thies, C. Landscape perspectives on agricultural intensification and biodiversity—ecosystem service management. Ecol. Lett. 8, 857–874. https://doi.org/10.1111/j.1461-0248.2005.00782.x (2005).
    Article  Google Scholar 

    39.
    Kiviniemi, K. & Eriksson, O. Size-related deterioration of semi-natural grassland fragments in Sweden. Divers. Distrib. 8, 21–29 (2002).
    Article  Google Scholar 

    40.
    Cousins, S. A. O. & Lindborg, R. Remnant grassland habitats as source communities for plant diversification in agricultural landscapes. Biol. Conserv. 141, 233–240. https://doi.org/10.1016/j.biocon.2007.09.016 (2008).
    Article  Google Scholar 

    41.
    Knapp, M. & Rezac, M. Even the smallest non-crop habitat islands could be beneficial: distribution of carabid beetles and spiders in agricultural landscape. PLoS ONE 10, e0123052. https://doi.org/10.1371/journal.pone.0123052 (2015).
    CAS  Article  PubMed  PubMed Central  Google Scholar 

    42.
    Oksuz, D. P. et al. Increasing biodiversity in wood-pastures by protecting small shrubby patches. For. Ecol. Manag. 464, 118041. https://doi.org/10.1016/j.foreco.2020.118041 (2020).
    Article  Google Scholar 

    43.
    Herzon, I. & Mikk, M. Farmers’ perceptions of biodiversity and their willingness to enhance it through agri-environment schemes: a comparative study from Estonia and Finland. J. Nat. Conserv. 15, 10–25 (2007).
    Article  Google Scholar 

    44.
    Wu, P. et al. Contrasting effects of natural shrubland and plantation forests on bee assemblages at neighboring apple orchards in Beijing, China. Biol. Conserv. 237, 456–462. https://doi.org/10.1016/j.biocon.2019.07.029 (2019).
    Article  Google Scholar 

    45.
    Liu, Y., Axmacher, J. C., Wang, C., Li, L. & Yu, Z. Ground beetles (Coleoptera: Carabidae) in the intensively cultivated agricultural landscape of Northern China—implications for biodiversity conservation. Insect Conserv. Divers. 3, 34–43. https://doi.org/10.1111/j.1752-4598.2009.00069.x (2010).
    Article  Google Scholar 

    46.
    Hodge, I. & Reader, M. The introduction of Entry Level Stewardship in England: extension or dilution in agri-environment policy?. Land Use Policy 27, 270–282. https://doi.org/10.1016/j.landusepol.2009.03.005 (2010).
    Article  Google Scholar 

    47.
    Landis, D. A. Designing agricultural landscapes for biodiversity-based ecosystem services. Basic Appl. Ecol. 18, 1–12. https://doi.org/10.1016/j.baae.2016.07.005 (2017).
    Article  Google Scholar 

    48.
    Tscharntke, T. et al. Landscape moderation of biodiversity patterns and processes—eight hypotheses. Biol. Rev. Camb. Philos. Soc. 87, 661–685. https://doi.org/10.1111/j.1469-185X.2011.00216.x (2012).
    Article  PubMed  Google Scholar 

    49.
    ESRI. ArcGIS 10.2 for Desktop. (2014).

    50.
    Han, Y. Study on Spatial and Temporal Patterns of Biodiversity in Intensive Agricultural Landscape of Quzhou County Master thesis (China Agricultural University, 2015).

    51.
    Hurlbert, A. H. Species-energy relationships and habitat complexity in bird communities. Ecol. Lett. 7, 714–720. https://doi.org/10.1111/j.1461-0248.2004.00630.x (2004).
    Article  Google Scholar 

    52.
    Zou, Y. et al. Diversity patterns of ground beetles and understory vegetation in mature, secondary, and plantation forest regions of temperate northern China. Ecol. Evol. 5, 531–542. https://doi.org/10.1002/ece3.1367 (2015).
    Article  PubMed  PubMed Central  Google Scholar 

    53.
    Jari, O. et al. In Package‘vegan’, Community Ecology Package 221 (2019).

    54.
    Faith, D. P., Minchin, P. R. & Belbin, L. Compositional dissimilarity as a robust measure of ecological distance. Vegetatio 69, 57–68. https://doi.org/10.1007/bf00038687 (1987).
    Article  Google Scholar 

    55.
    Ripley, B. et al. In Support Functions and Datasets for Venables and Ripley’s MASS (2019).

    56.
    Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A. & Smith, G. M. Mixed Effects Models and Extensions in Ecology with R. (2009).

    57.
    Roger, B. et al. In Spatial Dependence: Weighting Schemes, Statistics and Models 131–133 (2018).

    58.
    W. N. Venables, D. M. Smith & the_R_Core_Team. In Notes on R: A Programming Environment for Data Analysis and Graphics (Version 4.0.2, 2020).

    59.
    Lin, S. Honey Source Plant (China Forestry Publishing House, 1989).

    60.
    Xu, W. Chinese Honey Source Plant (Heilongjiang Science and Technology Press, 1983).

    61.
    Ke, X. Honey Powder Botanical (China Agriculture Press, 1995).

    62.
    Ma, D. & Liang, S. Chinese Honey Powder Source Plants and Their Utilization (1993). More

  • in

    On the robustness of an eastern boundary upwelling ecosystem exposed to multiple stressors

    1.
    Chavez, F. P. & Messié, M. A comparison of eastern boundary upwelling ecosystems. Prog. Oceanogr. 83, 80–96 (2009).
    ADS  Article  Google Scholar 
    2.
    Auger, P.-A., Gorgues, T., Machu, E., Aumont, O. & Brehmer, P. What drives the spatial variability of primary productivity and matter fluxes in the north-west African upwelling system? A modelling approach. Biogeosciences 13, 6419–6440 (2016).
    ADS  Article  Google Scholar 

    3.
    Benazzouz, A. et al. An improved coastal upwelling index from sea surface temperature using satellite-based approach—The case of the Canary Current upwelling system. Cont. Shelf Res. 81, 38–54 (2014).
    ADS  Article  Google Scholar 

    4.
    Citeau, J., Finaud, L., Cammas, J. & Demarcq, H. Questions relative to ITCZ migrations over the tropical Atlantic ocean, sea surface temperature and Senegal River runoff. Meteorol. Atmos. Phys. 41, 181–190 (1989).
    ADS  Article  Google Scholar 

    5.
    Maloney, E. & Shaman, J. Intraseasonal variability of the West African Monsoon and Atlantic ITCZ. J. Clim. 21, 2898–2918 (2008).
    ADS  Article  Google Scholar 

    6.
    Herbland, A. & Voituriez, B. L. production primaire dans l’upwelling mauritanien en mars 1973. Cahiers ORSTOM 12, 187–201 (1974).
    CAS  Google Scholar 

    7.
    Poloczanska, E. S. et al. Responses of marine organisms to climate change across oceans. Front. Mar. Sci. 3, 62 (2016).
    Article  Google Scholar 

    8.
    Hjermann, D. Ø., Ottersen, G. & Stenseth, N. C. Competition among fishermen and fish causes the collapse of Barents Sea capelin. PNAS 101, 11679–11684 (2004).
    ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

    9.
    Fréon, P., Cury, P., Shannon, L. & Roy, C. Sustainable exploitation of small pelagic fish stocks challenged by environmental and ecosystem changes: A review. Bull. Mar. Sci. 76, 385–462 (2005).
    Google Scholar 

    10.
    Schwartzlose, R. A. et al. Worldwide large-scale fluctuations of sardine and anchovy populations. Afr. J. Mar. Sci. 21, 289–347 (1999).
    Article  Google Scholar 

    11.
    Hofstede, R. T., Dickey-Collas, M., Mantingh, I. T. & Wague, A. The link between migration, the reproductive cycle and condition of Sardinella aurita off Mauritania, north-west Africa. J. Fish Biol. 71, 1293–1302 (2007).
    Article  Google Scholar 

    12.
    Zeeberg, J., Corten, A., Tjoe-Awie, P., Coca, J. & Hamady, B. Climate modulates the effects of Sardinella aurita fisheries off Northwest Africa. Fish. Res. 1, 65–75 (2008).
    Article  Google Scholar 

    13.
    Gibson, R. N., Atkinson, R. J. A. & Gordon, J. D. M. Oceanography and Marine Biology, Vol. 47 (2009).

    14.
    Hays, G. C., Richardson, A. J. & Robinson, C. Climate change and marine plankton. Trends Ecol. Evol. 20, 337–344 (2005).
    PubMed  Article  PubMed Central  Google Scholar 

    15.
    Ndoye, S. et al. Dynamics of a “low-enrichment high-retention” upwelling center over the southern Senegal shelf. Geophys. Res. Lett. 44, 5034–5043 (2017).
    ADS  Article  Google Scholar 

    16.
    Behagle, N. et al. Acoustic distribution of discriminated micronektonic organisms from a bi-frequency processing: The case study of eastern Kerguelen oceanic waters. Prog. Oceanogr. 156, 276–289 (2017).
    Article  Google Scholar 

    17.
    Benoit-Bird, K. & Au, W. Diel migration dynamics of an island-associated sound-scattering layer. Deep Sea Res. Part I 51, 707–719 (2004).
    Article  Google Scholar 

    18.
    Sato, M. & Benoit-Bird, K. J. Spatial variability of deep scattering layers shapes the Bahamian mesopelagic ecosystem. Mar. Ecol. Prog. Ser. 580, 69–82 (2017).
    ADS  Article  Google Scholar 

    19.
    Algueró-Muñiz, M. et al. Ocean acidification effects on mesozooplankton community development: Results from a long-term mesocosm experiment. PLoS ONE 12, e0175851 (2017).
    PubMed  PubMed Central  Article  CAS  Google Scholar 

    20.
    Matlab R 2018a. The Math Works, Inc. (MATLAB & Simulink – MathWorks., 2018).

    21.
    Hegerl, G. C. et al. Causes of climate change over the historical record. Environ. Res. Lett. 14, 123006 (2019).
    ADS  CAS  Article  Google Scholar 

    22.
    Belkin, I. M. Rapid warming of large marine ecosystems. Prog. Oceanogr. 81, 207–213 (2009).
    ADS  Article  Google Scholar 

    23.
    Valdés, L. & Déniz-González, I. Oceanographic and biological features in the Canary Current Large Marine Ecosystem, Vol. 115 (2015).

    24.
    Gómez-Letona, M., Ramos, A. G., Coca, J. & Arístegui, J. Trends in primary production in the canary current upwelling system—A regional perspective comparing remote sensing models. Front. Mar. Sci. 4, 370 (2017).
    Article  Google Scholar 

    25.
    Bakun, A. Global climate change and intensification of coastal ocean upwelling. Science 247, 198–201 (1990).
    ADS  CAS  PubMed  Article  Google Scholar 

    26.
    Benazzouz, A., Demarcq, H. & González-Nuevo, G. Oceanographic and biological features in the Canary current large marine ecosystem. (2015).

    27.
    Behrenfeld, M. J. et al. Climate-driven trends in contemporary ocean productivity. Nature 444, 752–755 (2006).
    ADS  CAS  PubMed  Article  Google Scholar 

    28.
    Boyce, D. G., Lewis, M. R. & Worm, B. Global phytoplankton decline over the past century. Nature 466, 591–596 (2010).
    ADS  CAS  PubMed  Article  Google Scholar 

    29.
    Hofmann, M., Worm, B., Rahmstorf, S. & Schellnhuber, H. J. Declining ocean chlorophyll under unabated anthropogenic CO2 emissions. Environ. Res. Lett. 6, 034035 (2011).
    ADS  Article  CAS  Google Scholar 

    30.
    Lewandowska, A. et al. Effects of sea surface warming on marine plankton. Ecol. Lett. 17, 614–623 (2014).
    PubMed  Article  PubMed Central  Google Scholar 

    31.
    Arístegui, J. et al. Sub-regional ecosystem variability in the Canary Current upwelling. Prog. Oceanogr. 83, 33–48 (2009).
    ADS  Article  Google Scholar 

    32.
    Demarcq, H. Trends in primary production, sea surface temperature and wind in upwelling systems (1998–2007). Prog. Oceanogr. 83, 376–385 (2009).
    ADS  Article  Google Scholar 

    33.
    Barton, A. D., Irwin, A. J., Finkel, Z. V. & Stock, C. A. Anthropogenic climate change drives shift and shuffle in North Atlantic phytoplankton communities. PNAS 113, 2964–2969 (2016).
    ADS  CAS  PubMed  Article  Google Scholar 

    34.
    Jacob, B. G. et al. Major changes in diatom abundance, productivity, and net community metabolism in a windier and dryer coastal climate in the southern Humboldt Current. Prog. Oceanogr. 168, 196–209 (2018).
    ADS  Article  Google Scholar 

    35.
    Jacox, M. G., Hazen, E. L. & Bograd, S. J. optimal environmental conditions and anomalous ecosystem responses: Constraining bottom-up controls of phytoplankton biomass in the California current system. Sci. Rep. 6, 27612 (2016).
    ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

    36.
    Botsford, L. W., Lawrence, C. A., Dever, E. P., Hastings, A. & Largier, J. Wind strength and biological productivity in upwelling systems: An idealized study. Fish. Oceanogr. 12, 245–259 (2003).
    Article  Google Scholar 

    37.
    García-Reyes, M., Largier, J. L. & Sydeman, W. J. Synoptic-scale upwelling indices and predictions of phyto- and zooplankton populations. Prog. Oceanogr. 120, 177–188 (2014).
    ADS  Article  Google Scholar 

    38.
    Libralato, S., Coll, M., Tudela, S., Palomera, I. & Pranovi, F. Novel index for quantification of ecosystem effects of fishing as removal of secondary production. Mar. Ecol. Prog. Ser. https://doi.org/10.3354/meps07224 (2008).
    Article  Google Scholar 

    39.
    Gasol, J. M., del Giorgio, P. A. & Duarte, C. M. Biomass distribution in marine planktonic communities. Limnol. Oceanogr. 42, 1353–1363 (1997).
    ADS  CAS  Article  Google Scholar 

    40.
    Harfoot, M. B. J. et al. Emergent global patterns of ecosystem structure and function from a mechanistic general ecosystem model. PLoS Biol. 12, e1001841 (2014).
    PubMed  PubMed Central  Article  CAS  Google Scholar 

    41.
    Wang, H., Morrison, W., Singh, A. & Weiss, H. General Mechanisms for Inverted Biomass Pyramids in Ecosystems. arXiv:0811.3657. [q-bio] (2008).

    42.
    Benoit-Bird, K. J. & Lawson, G. L. Ecological insights from pelagic habitats acquired using active acoustic techniques. Ann. Rev. Mar. Sci. 8, 463–490 (2016).
    PubMed  Article  PubMed Central  Google Scholar 

    43.
    Alcaraz, M., Felipe, J., Grote, U., Arashkevich, E. & Nikishina, A. Life in a warming ocean: Thermal thresholds and metabolic balance of arctic zooplankton. J. Plankton Res. 36, 3–10 (2014).
    Article  Google Scholar 

    44.
    Brochier, T. et al. Complex small pelagic fish population patterns arising from individual behavioral responses to their environment. Prog. Oceanogr. 164, 12–27 (2018).
    ADS  Article  Google Scholar 

    45.
    Richardson, A., Schoeman, D., Richardson, A. J. & Schoeman, D. S. Climate impact on plankton ecosystems in the Northeast Atlantic. Science 305, 1609–1612 (2004).
    ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

    46.
    Braham, C.-B. & Corten, A. Pelagic fish stocks and their response to fisheries and environmental variation in the Canary Current large marine ecosystem. Oceanographic and biological features in the Canary Current Large Marine Ecosystem 197–213 (2015).

    47.
    Ba, K. et al. Resilience of key biological parameters of the Senegalese flat sardinella to overfishing and climate change. PLoS ONE 11, e0156143 (2016).
    PubMed  PubMed Central  Article  CAS  Google Scholar 

    48.
    Thiaw, M. et al. Effect of environmental conditions on the seasonal and inter-annual variability of small pelagic fish abundance off North-West Africa: The case of both Senegalese sardinella. Fish. Oceanogr. https://doi.org/10.1111/fog.12218 (2017).
    Article  Google Scholar 

    49.
    Sarré, A. et al. Climate-driven shift of Sardinella aurita stock in Northwest Africa ecosystem as evidenced by robust spatial indicators [résumé]. In International conference ICAWA 2016 : extended book of abstract : the AWA project : ecosystem approach to the management of fisheries and the marine environment in West African waters (eds. Brehmer, P. et al.) 67–68 (SRFC/CSRP, 2017).

    50.
    Richardson, A. J. In hot water: Zooplankton and climate change. ICES J. Mar. Sci. 65, 279–295 (2008).
    Article  Google Scholar 

    51.
    Beaugrand, G., Reid, P., Ibañez, F., Lindley, J. & Edwards, M. Reorganization of North Atlantic marine copepod biodiversity and climate. Science 296, 1692–1694 (2002).
    ADS  CAS  PubMed  Article  Google Scholar 

    52.
    Lindley, J. A. & Daykin, S. Variations in the distributions of Centropages chierchiae and Temora stylifera (Copepoda: Calanoida) in the north-eastern Atlantic Ocean and western European shelf waters. ICES J. Mar. Sci. 62, 869–877 (2005).
    Article  Google Scholar 

    53.
    Chassot, E. et al. Global marine primary production constrains fisheries catches. Ecol. Lett. 13, 495–505 (2010).
    PubMed  Article  Google Scholar 

    54.
    Diogoul, N. et al. Fine-scale vertical structure of sound-scattering layers over an east border upwelling system and its relationship to pelagic habitat characteristics. Ocean Sci. 16, 65–81 (2020).
    ADS  CAS  Article  Google Scholar 

    55.
    Brehmer, P. et al. Schooling behaviour of small pelagic fish: Phenotypic expression of independent stimuli. Mar. Ecol. Prog. Ser. 334, 263–272 (2007).
    ADS  Article  Google Scholar 

    56.
    Munday, P. L., Jones, G. P., Pratchett, M. S. & Williams, A. J. Climate change and the future for coral reef fishes. Fish Fish. 9, 261–285 (2008).
    Article  Google Scholar 

    57.
    Costello, J. H., Sullivan, B. K. & Gifford, D. J. A physical–biological interaction underlying variable phenological responses to climate change by coastal zooplankton. J. Plankton Res. 28, 1099–1105 (2006).
    Article  Google Scholar 

    58.
    Garzke, J., Ismar, S. M. H. & Sommer, U. Climate change affects low trophic level marine consumers: Warming decreases copepod size and abundance. Oecologia 177, 849–860 (2015).
    ADS  PubMed  Article  PubMed Central  Google Scholar 

    59.
    Horne, C. R., Hirst, A. G., Atkinson, D., Neves, A. & Kiørboe, T. A global synthesis of seasonal temperature–size responses in copepods. Glob. Ecol. Biogeogr. 25, 988–999 (2016).
    Article  Google Scholar 

    60.
    Clark, C. W. & Levy, D. A. Diel vertical migrations by juvenile sockeye salmon and the antipredation window. Am. Nat. 131, 271–290 (1988).
    Article  Google Scholar 

    61.
    Lampert, W. The adaptive significance of diel vertical migration of zooplankton. Funct. Ecol. 3, 21–27 (1989).
    Article  Google Scholar 

    62.
    Hansson, S. Variation in hydroacoustic abundance of pelagic fish. Fish. Res. 16, 203–222 (1993).
    Article  Google Scholar 

    63.
    Tiedemann, M. & Brehmer, P. Larval fish assemblages across an upwelling front: Indication for active and passive retention. Estuar. Coast. Shelf Sci. 187, 118–133 (2017).
    ADS  Article  Google Scholar 

    64.
    CCLME. Analyse diagnostique transfrontalière du Grand écosystème marin du courant des Canaries 144 (2016).

    65.
    Bernhardt, J. R. & Leslie, H. M. Resilience to climate change in coastal marine ecosystems. Annu. Rev. Mar. Sci. 5, 371–392 (2013).
    Article  Google Scholar 

    66.
    Baldé, B. S. et al. Variability of key biological parameters of round sardinella Sardinella aurita and the effects of environmental changes. J. Fish Biol. 94, 391–401 (2019).
    PubMed  Article  PubMed Central  Google Scholar 

    67.
    Binet, D. Rôle possible d’une intensification des alizés sur le changement de répartition des sardines et sardinelles le long de la côte Ouest africaine. Aquat. Living Resour. 1, 115–132 (1988).
    Article  Google Scholar 

    68.
    Berraho, A., Somoue, L., Hernández‐León, S. & Valdés, L. Zooplankton in the canary current large marine ecosystem. In Oceanographic and biological features in the Canary Current Large Marine Ecosystem Vol. 115, 183‐195 (IOC Technical Series, 2015).

    69.
    Ndour, I., Berraho, A., Fall, M., Ettahiri, O. & Sambe, B. Composition, distribution and abundance of zooplankton and ichthyoplankton along the Senegal-Guinea maritime zone (West Africa). Egypt. J. Aquat. Res. 44, 109–124 (2018).
    Article  Google Scholar 

    70.
    Sarré, A., Krakstad, J.-O., Brehmer, P. & Mbye, E. M. Spatial distribution of main clupeid species in relation to acoustic assessment surveys in the continental shelves of Senegal and The Gambia. Aquat. Living Resour. 31, 9 (2018).
    Article  Google Scholar 

    71.
    Foote, K. G., Knudsen, H. P., Vestnes, G., MacLennan, D. N. & Simmonds, E. J. Technical Report: ‘“Calibration of acoustic instruments for fish density estimation: A practical guide”’. J. Acoust. Soc. Am. 83, 831–832 (1987).
    Google Scholar 

    72.
    Perrot, Y. et al. Matecho: An open-source tool for processing fisheries acoustics data. Acoust. Aust. 46, 241–248 (2018).
    Article  Google Scholar 

    73.
    MacLennan, D. N., Fernandes, P. G. & Dalen, J. A consistent approach to definitions and symbols in fisheries acoustics. ICES J. Mar. Sci. 59, 365–369 (2002).
    Article  Google Scholar 

    74.
    Jech, J. M., Lawson, G. L. & Lavery, A. C. Wideband (15–260 kHz) acoustic volume backscattering spectra of Northern krill (Meganyctiphanes norvegica) and butterfish (Peprilus triacanthus). ICES J. Mar. Sci. 74, 2249–2261 (2017).
    Article  Google Scholar 

    75.
    Madureira, L. S. P., Everson, I. & Murphy, E. J. Interpretation of acoustic data at two frequencies to discriminate between Antarctic krill (Euphausia superba Dana) and other scatterers. J. Plankton Res. 15, 787–802 (1993).
    Article  Google Scholar 

    76.
    Brehmer, P., Georgakarakos, S., Josse, E., Trygonis, V. & Dalen, J. Adaptation of fisheries sonar for monitoring schools of large pelagic fish: Dependence of schooling behaviour on fish finding efficiency. Aquat. Living Resour. 20, 377–384 (2007).
    Article  Google Scholar 

    77.
    D’Elia, L. et al. A longitudinal study of the teacch program in different settings: The potential benefits of low intensity intervention in preschool children with autism spectrum disorder. J. Autism Dev. Disord. 44, 615–626 (2014).
    PubMed  Article  Google Scholar 

    78.
    Zwolinski, J., Morais, A., Marques, V., Stratoudakis, Y. & Fernandes, P. G. Diel variation in the vertical distribution and schooling behaviour of sardine (Sardina pilchardus) off Portugal. ICES J. Mar. Sci. 64, 963–972 (2007).
    Article  Google Scholar 

    79.
    Ayoubi, S. E. et al. Estimation of target strength of Sardina pilchardus and Sardinella aurita by theoretical approach. Fish. Sci. 82, 417–423 (2016).
    Article  CAS  Google Scholar 

    80.
    Saunders, R. A., Fielding, S., Thorpe, S. E. & Tarling, G. A. School characteristics of mesopelagic fish at South Georgia. Deep Sea Res. Part I 81, 62–77 (2013).
    Article  Google Scholar 

    81.
    R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2016).

    82.
    Suppiah, R. & Hennessy, K. J. Trends in total rainfall, heavy rain events and number of dry days in Australia, 1910–1990. Int. J. Climatol. 18, 1141–1164 (1998).
    Article  Google Scholar 

    83.
    Cotter, J. A selection of nonparametric statistical methods for assessing trends in trawl survey indicators as part of an ecosystem approach to fisheries management (EAFM). Aquat. Living Resour. 22, 173–185 (2009).
    Article  Google Scholar  More

  • in

    Biting midge dynamics and bluetongue transmission: a multiscale model linking catch data with climate and disease outbreaks

    We have used a multi-scale modelling approach combining multiple modelling and inference types and techniques: generalised linear mixed-effect models (GLMMs), a mechanistic BTV transmission simulation model, marginal maximum likelihood inference, and direct calculation of resultant BTV risk predictions. We have parameterised our models with data from multiple existing sources (spatio-temporal climate time series, livestock density estimates, midge life process rates, BTV incidence time series), as well as data from our recent field study on the activity of midges at different latitudes in Europe34. The workflow and data sources of this study are summarised in Fig. 1.
    Biting midge collection and identification
    The method of biting midge sampling and identification was described in our earlier study34. In short, three habitat types were defined in which Culicoides specimens were collected: “farm”, “peri-urban”, and “wetland”. Traps were placed within a 50 m radius of cattle, a house, or waterbody, for farm, peri-urban, and wetland habitats respectively. Habitat types generally matched the classification of the CORINE European Land cover database46. Collections were performed in Sweden (surroundings of Linköping N58.410808, E15.621532), The Netherlands (surroundings of Wageningen N51.964795, E5.662898), and Italy (surroundings of San Benedetto del Tronto N42.949483, E13.878503). Onderstepoort Veterinary Institute black light traps were placed at three independent locations for each selected habitat type. Traps were at least 100 m apart to prevent overlap of the active trapping radius47. More details and exact trap locations can be found in Vogels et al.48 and Möhlmann et al.34. Collections were performed for six consecutive days in each month in all three countries, during the period from July 2014 to June 2015 except the winter months of December, January and February (and March for Sweden). Traps were emptied and repositioned at different locations every 24 h. Collections were sorted and stored in 70% ethanol at − 20 °C. Samples were identified to species level using the Interactive Identification Key for Culicoides (IIKC) developed by Mathieu et al.49.
    Predictor variables for midge activity
    Tinytag® meteorological data loggers (Gemini Data Loggers, Chichester, UK) were used to record local temperature and relative humidity every 30 min during the collection period from 17th July 2014 until 3rd July 2015. For each habitat in each country, one data logger was used (3 countries × 3 habitats). Additional meteorological data was collected from a number of sources: hourly wind velocities and local temperatures were available from the weather station closest to the trap location in Sweden (Swedish Meteorological and Hydrological Institute (SMHI) Linköping weather station), The Netherlands (Koninklijk Nederlands Meteorologisch Instituut (KNMI) weather stations “De Bilt”, “Deelen”, “Cabauw”, and “Volkel”), and Italy (San Benedetto del Tronto weather station). Finally, additional minimum, maximum, and mean daily temperatures along with precipitation and air pressure were sourced from the E-OBS European climate database50. This data was available at a daily temporal resolution and a spatial resolution of 0.25 degrees latitude and longitude. For climatic variables such as wind only one source per area could be used, whereas for temperature we had the opportunity to explore the most predictive of a number of data sources per country. Habitat effect was included by using the habitat types (farm, peri-urban, wetland) as a categorical predictor in the regression models.
    BTV-seroprevalence survey data from Dutch sentinel herds
    After the initial BTV outbreaks in 2006, the Dutch government decided to establish a sentinel network of dairy cattle herds in the winter of 2006/2007 to monitor the re-emergence of BTV-8 in 2007, using repeated milk ELISA testing35,36. For study purposes, The Netherlands was divided into 20 compartments based on geographic boundaries as proposed in Commission Decision 2005/393/EC. In each compartment, at least 10 randomly selected herds had to be sampled (with at least sixteen cows per herd) to obtain the required sample size. Herds were not necessarily completely BTV-8 seronegative at initial investigation, but cows designated for the sentinel program had to be BTV-8 seronegative at the moment of selection in May 2007. Therefore, dairy herds to act as sentinels for BTV incidence were selected, that had at least sixteen seronegative cows and at least 50 cattle in total. Monthly milk samples were collected from the sentinel cows in each herd unless prevalence had already reached 100%. The first round of monthly testing of sentinel cows was done in June 2007 and continued until January 2008. The monthly milk samples were tested at GD Animal Health, Deventer The Netherlands for antibodies to BTV-8 using a commercially available ELISA test. For further details on the sampling protocol and commercial ELISA see Santman-Berends et al.35,36.
    Regression model for biting midge catches
    A preliminary inference investigation using a hurdle negative binomial model to explain trap catches, inferred using Metropolis–Hastings MCMC, suggested that the excess of zero catches could be explained without zero inflation using the hurdle mechanism. We therefore used GLMM regression for its convenience and flexibility. Regression for the fixed effect coefficients and variance parameters of the random effects was performed via maximum likelihood using the Laplace approximation method implemented by the fitglme MATLAB® function. We followed the recommendation guidelines of Bolker et al.51 for using generalized linear models in the context of applied ecology, starting from a model with a full set of predictors and performed systematic model reduction using AICc to score model improvement (backwards model selection). AICc is a well-established information criterion for model selection since it is easy to calculate and interpret for GLMMs. The AICc difference between two models (ΔAICc) estimates the relative likelihood of the two models (with (sim {text{exp(}} – Delta {text{AICc}}/2)) as relative likelihood for the model with greater AICc). Moreover, the model selected by lowest AICc is also the model that would be selected by leave-one-out cross-validation52.
    Candidate predictor variables for removal from the model were chosen by assessing which fixed-effect coefficients had the greatest P value (for the null hypothesis that the coefficient is zero) and which random effect coefficients had the smallest predicted standard deviation. This was followed by trialling the removal of either of the variables and removing the variable that reduced AICc the most, until no further reduction could be made. Several trials were made where we started from a number of different highly over-parameterised models, which all ended with the same best model.
    For any given combination of predictor variables, the catches were assumed to be conditionally Poisson distributed, with the conditional mean for each collection defined by the log-link relationship,

    $$lnleft( {mu_{lct} } right) = beta cdot X_{lt} + b_{l} cdot Z_{lt} + rho_{ct} + varepsilon_{lt}$$
    (1)

    where (beta) is the fixed effect regression coefficient that applies to all catches, with ({X}_{lt}) denoting the fixed effect predictors at trap location l on day t. ({b}_{l}sim Normal(0,Sigma )) are the location-grouped random effect regression coefficients with covariance matrix (Sigma), with ({Z}_{lt}) denoting the random effect predictors at location l on day t. The number of midges per catch was highly variable. Therefore, we included an independent random effect for each catch location and day ({varepsilon }_{lt}sim Normal(0,{sigma }_{varepsilon })). This corresponds to assuming that the number of midges per catch is Poisson log-normally distributed (a standard distribution for over-dispersed count data53). Spatial autocorrelation has been found in other midge catch regression studies38, so we also included an autocorrelation random effect grouped by country of trapping and day; ({rho }_{ct}sim Normal(0,{sigma }_{varrho })). This effect accounts for events influencing trap catch at a wider scale than just the location definition of each trap. See supplementary information for full model variables, regression coefficients and random effect variance estimates as well as AICc improvements from other models.
    Connecting the regression model for midge captures to a midge biting prediction for herds
    We connected the regression model for daily midge trap catch to a prediction of daily biting on cattle. In the BTV modelling literature it is common to assume that the vector-to-host ratio, and therefore the biting rate per livestock host, is proportional to the expected capture rate of a midge trap catch regression model28,54. In this study, we assumed that the biting rate per livestock host is proportional to the regression model presented in Results, with the major difference that we infer the scaling parameter that best fitted the serological data from Dutch sentinel herds. The location-group random effects in the regression model modelled unexplained spatial variation in average midge capture counts between trapping locations. We assumed that this unexplained spatial variation in midge abundance (as measured by trapping) mirrored the spatial variation in midge biting between different herds. Combining the proportionality assumption with the spatial variation assumption gave the biting rate among herds as,

    $${text{Biting }},{text{rate }},{text{of }},{text{susceptible}},{text{ midges }},{text{per}},{text{ cattle }},{text{at }},{text{herd}},h,{text{on}},{text{ day}},t = xi mu_{ht} .$$
    (2)

    where ({mu }_{ht}) is the expected trap catch from the regression model given the climate condition local to herd h on day t and (xi) is a scaling parameter between the mean catch prediction and the biting rate prediction. In order to infer the proportionality factor between the catch model and the daily biting rate ((xi)), the outcomes of a dynamic transmission model for each herd were linked to the observed BTV seroprevalence data (see—“Mechanistic transmission model for BTV transmission within herds” section).
    Mechanistic transmission model for BTV transmission within herds
    We used a dynamic and mechanistic model of BTV transmission within herds, which was then matched to the data from the Dutch sentinel study. The dynamic BTV transmission model was formulated using disease compartments and rate-based transitions (see Keeling and Rohani55 for further details on this class of transmission model). In addition, it is in most respects similar to the model presented by Gubbins et al.27 in treating infectiousness amongst cattle, and latency amongst midges, as multi-stage processes that evolve deterministically. In particular, we follow Gubbins et al. in modelling the life processes of the infectious and latent midges (mortality, extrinsic incubation of BTV, and biting) as temperature dependent, and therefore varying daily with local background atmospheric temperature. (see Fig. 6 for a schematic diagram of the transmission model).
    Figure 6

    Schematic representation of the cattle herd level BTV transmission model. The population of cattle and infected biting midges are divided amongst discrete disease compartments. BTV latent midges (EM) enter the model at a rate proportional to the daily prediction of the catch model. The extrinsic incubation period for latent midges is modelled as a multi-stage process before midges become infectious (IM). Susceptible cattle (SC) become infectious cattle (IC) after a bite from infectious midges (IM), to become resistant cattle (RC) in time (also modelled as a multi-stage process; red box). Transitions are shown as solid lines, coloured according to their dependence on environmental variables: constant per-capita (black), daily mean temperature dependent (red), all predictor variables of capture model and the catch-to-bite scale parameter ({varvec{xi}}) (blue). Dotted lines indicate where the number of infected individuals in one species increases the incidence rate in the other species. Outcomes of the model are linked to observed cattle milk serology time series by the first two infectious stages for cattle with virus being undetectable by ELISA (blue box), whereas subsequent infectious stages and the recovered stage are detectable by ELISA (green box). The likelihood function for (xi) was inferred by marginalisation over the latent stochastic variables affecting model outcomes (e.g. herd-specific random effects, daily fluctuations in midge activity). Used images were available under open licence Creative Commons Deed CC0.

    Full size image

    The dynamic and mechanistic BTV transmission model used in this study describes the evolution of the numbers of susceptible, infected, and recovered cattle as well as latent and infectious midges for each herd (Fig. 6). Temperature-dependent midge bionomic rates were used for biting frequency of individual infectious midges56, the incubation rate of BTV within the midge57, and for the midge mortality58. The midge bionomic rates at each herd on each day were determined by the local mean temperature day according to the E-OBS climate dataset (see above—Predictor variables for midge activity). We modelled the incubation period of BTV within the midge vector as a ten-stage process, which is within the range of best fit models found in meta-analysis and laboratory studies of BTV incubation57 (see supplementary information for complete model details and literature estimates for rates).
    The period during which BTV-infected cattle are detectable using an ELISA test (typically from 8–9 days post-infection onwards59) does not match the period during which the cattle are infectious (rapidly post-infection and then for an average of 20.6 days27). BTV-infected cattle can be in four states that are relevant to transmission modelling and their milk serology: 1) uninfected and susceptible to BTV, 2) infectious but undetectable by milk ELISA, 3) infectious and detectable by milk ELISA or, 4) non-infectious recovered from BTV but still milk ELISA positive. The BTV infectious period for cattle is usually modelled as a 5-stage process27, therefore it was convenient to model cattle in the first two stages of their infectious period (an average duration of 8.2 days) as infectious but undetectable. Cattle in the final three stages of the BTV infectious period are infectious and detectable (see Fig. 6 for a schematic representation of the BTV transmission and serology model).
    The rate at which new susceptible midges arrived to bite cattle was assumed to be proportional to the expected number of midges from the regression model (({mu }_{ht})), with predictions of midge catches on each day and in each regional compartment of The Netherlands whilst the sentinel herd study was ongoing using the E-OBS historic climate records (see above—“Connecting the regression model for midge captures to a midge biting prediction for herds” section). The random effects in the catch model imply that ({mu }_{ht}) is a daily varying random variable, and that our transmission model is in the class of piecewise-deterministic Markov processes60. We assumed that the location-grouped random effects observed in the catch model became herd-grouped random effects for the biting model. In other words, we assumed that the high variance in midge catching between trapping location reflects high variance in midge biting between different cattle herd locations. Although an assumption, this would explain the highly variable intensity of BTV transmission observed between different herds in The Netherlands sentinel survey despite each herd experiencing a similar climate36. Because the herd locations were known only as geographic compartment occupancy, the daily local climate variables from the gridded E-OBS data used for predicting ({mu }_{ht}) were averaged over all spatial grids overlapping the herd’s geographic compartment. Accessing detailed and spatio-temporally resolved wind data across Europe was challenging, therefore we used the long-term average wind velocity of The Netherlands weather stations (see above—“Predictor variables for biting midge activity” section) as a constant predictor.
    When simulating an outbreak of BTV within a herd, we first determined all relevant climatic predictors for the herd’s regional compartment and the daily temperature dependent midge bionomic rates. Second, we generated the herd-grouped and daily varying random effect coefficients which determined how biting at the herd from susceptible midges varied from a median prediction. Third, we solved the resultant deterministic BTV transmission model for each farm using the ode45 MATLAB® function (see Fig. 6 for an overview).
    Inferring the catch-to-biting scale parameter from serological data
    We inferred a maximum likelihood estimator for the trap-to-bite scale parameter by repeated simulation of the percentage of cattle detectable by milk ELISA tests herds in the sentinel herd network. For this we used the climatic conditions of The Netherlands in 2007, and at each simulation repeated redrawing the unobserved random effects for each herd and day. The average likelihood over many repeated simulations corresponds to a Monte Carlo estimate of the true marginal likelihood of the parameter (xi). Estimating the likelihood over a range of values of (xi) allowed the construction of a log-likelihood profile.
    The stochastic elements of the piecewise-deterministic BTV transmission model were (i) the herd-grouped random coefficients (this modelled how biting varied from herd-to-herd) and (ii) the daily varying random effects (this modelled how biting varied from day-to-day). It is convenient to denote ({W}_{h}= ({{b}_{l}}^{(h)},{{rho }_{0}}^{(h)},{{rho }_{1}}^{(h)},{{rho }_{2}}^{(h)},…,{{epsilon }_{0}}^{(h)},{{epsilon }_{1}}^{(h)},{{epsilon }_{2}}^{(h)},…)) as the collection of all stochastic elements for herd h. For each simulation of the transmission model in each herd, we first drew ({W}_{h}) from their inferred distribution (see supplementary information for distribution parameters of best fitting model).
    The likelihood of ({W}_{h}) and (xi) for each herd h was the chance of selecting the numbers of ELISA seroconverted cattle observed at the herd each month by The Netherlands sentinel study from the underlying distribution of ELISA detectable cattle implied by simulating the transmission model conditional on ((xi ,{W}_{h})),

    $$L_{h} left( {xi ,W_{h} } right) = P({text{Serology }},{text{data }},{text{collected }},{text{at }},{text{herd }},h| xi ,W_{h} ).$$
    (3)

    Since we were not interested in inferring the specific values of ({W}_{h}) for each herd, they were treated as “nuisance” parameters. We inferred a maximum likelihood estimate, with confidence intervals, for (xi) by first estimating the marginal likelihood for (xi) (that is the likelihood after integrating over all possible values of the nuisance parameters) at each herd h,

    $${L}_{h}(xi ) = int {L}_{h}left(xi ,wright)fleft(wright)dw.$$
    (4)

    where (f) is the density function for the distribution of random effects derived from the trapping model. The marginal log-likelihood function for the trap-to-bite scaling parameter, (l(xi )), for serological data over a number of herds, was then just the sum of the individual herd marginal log-likelihoods,

    $$l(xi ) = sum_{h}log {L}_{h}left(xi right).$$
    (5)

    The herds we chose to contribute to the log-likelihood were those where BTV was found to be already present at the beginning of the study (see supplementary information for more details), to avoid making further assumptions about the introduction mechanism into the herd.
    In practice, the log-likelihood was estimated for a profile of values of (xi) by simulating multiple realisations of ({W}_{h}) for each herd, that is we estimated (4) by Monte Carlo integration for (3) over a range of values of (xi), and interpolating between points with polynomial regression. The maximum likelihood estimator, ({xi }^{*}), was the maximizer of the marginal log-likelihood function presented in the main text along with confidence intervals derived by a standard comparison to the ({chi }^{2}) distribution (see methods section of King et al.61 for a brief but comprehensive introduction to maximum likelihood estimation using log-likelihood profiles in the context of inference for dynamical systems).
    Calculating and mapping the herd reproductive ratio for bluetongue
    The reproductive ratio for BTV will differ from day to day and across space. This reflects seasonality and variation in both climatic trends, and the population density of midges and livestock hosts. We approached estimating the reproductive ratio for BTV in the spirit of the case reproductive ratio62 using a technique already developed for midges spreading BTV37. That is, we calculated the expected number of secondary cases amongst hosts due to a host initially infected on each day t in each grid cell x whilst taking into account how the conditions for BTV transmission at location x changed after time t, and using the maximum likelihood model for midge biting. The size of each grid cell was determined by the resolution of the relevant climate data. We used the 0.25 degrees longitude and latitude grid resolution of the E-OBS climate datasets to map reproductive ratio estimates for Europe across space and time. Cattle and sheep densities across Europe were drawn from the livestock Geo-Wiki dataset63. The livestock Geo-Wiki datasets were more finely resolved (0.0083 degrees grid resolution) than the E-OBS climate datasets. We calculated the reproductive ratio for each grid cell x on each day t at the resolution of the E-OBS datasets. The cattle and sheep densities for this coarser grained grid were the average densities over the Geo-Wiki cells contained within the coarser grained grid.
    The average number of secondary BTV cases amongst all hosts (cattle and sheep) given a host initially infected on day t and at grid cell x, is denoted ({R}^{(C)}(x,t)) for an initial infected cow and ({R}^{(S)}(x,t)) for an initial infected sheep. For both host species, the average number of secondary cases can be calculated by considering; how many days the host’s viraemia will last, the rate at which the host is bitten each day, the percentage of the biting midges that will become infected, how many of these biting midges are expected to survive their EIP to become actively infectious, and how many livestock will be successfully infected by those actively infectious midges. The methodology for combining these estimates using information about midge bionomic rates, EIP distribution, and the daily temperatures on each day after the initial host was infected has been developed by Brand and Keeling37.
    In this study, we adapted the Brand-method for calculating the reproductive ratio to two species, and used the catch-to-biting scalar derived from comparison between the mechanistic transmission model and the herd sentinel serological survey. The cross-transmission between host species depends on how midge bites are distributed between cattle and sheep. We estimated the proportion of midge bites on cattle at grid cell x, ({phi }^{(C)}(x)), given the availability of sheep using a common relative preference model, e.g. Szmaragd et al.64,

    $${phi }^{left(Cright)}left(xright)=frac{{N}^{left(Cright)}left(xright)}{{N}^{left(Cright)}left(xright)+ pi {N}^{left(Sright)}left(xright)}.$$
    (6)

    where ({N}^{(C)}(x)) and ({N}^{(S)}(x)) are, the local density of cattle and sheep at grid cell x. Parameter (pi) is a measure of the vector preference for sheep compared to cattle; (pi 1) preference for sheep. A relative biting study for sheep and cattle has revealed a preference for biting cattle30, from which we derived an estimate (pi =0.115) for use in this study. We combined ({R}^{(C)}(x,t)) and ({R}^{(S)}(x,t)) into a single reproductive ratio by calculating the leading eigenvalue of the next-generation matrix65,

    $$Rleft(x,tright)= sqrt{{phi }^{left(Cright)}left(xright){R}^{left(Cright)}left(x,tright)+left(1-{phi }^{left(Cright)}left(xright)right){R}^{left(Sright)}left(x,tright)}.$$
    (7)

    An attractive feature of using the reproductive ratio as a measure of transmission intensity is its uncomplicated relationship with the persistence of transmission; if (Rle 1) then the infectious pathogen cannot persist. However, we expect that the biting rate, and therefore the reproductive ratio, will vary from herd-to-herd. From Eqs. (1) and (5) we see that the rate of biting from the susceptible midge population at each herd on each day depends on the random coefficients, ({{b}_{l}}^{(h)}), and daily varying random effects,({{rho }_{ct}}^{(h)}) and ({{varepsilon }_{t}}^{(h)}),

    $$text{Biting rate of susceptible midges per cattle at herd }htext{ on day }tpropto expleft({{b}_{l}}^{left(hright)}cdot {Z}_{lt} +{{rho }_{ct}}^{left(hright)} + {{varepsilon }_{t}}^{left(hright)} right).$$
    (8)

    The daily varying random effects ((rho) and (epsilon)) are averaged over our estimates for ({R}^{(C)}) and ({R}^{(S)}) (this can be achieved analytically; see supplementary information for further details), and therefore our estimate of the reproductive ratio does not depend on daily fluctuations in midge activity. However, variation in the herd-grouped random coefficients indicated systematic differences in midge activity between herds that will not ‘average out’ over time. The distribution of ({{b}_{l}}^{(h)}) therefore implied a distribution of biting rates for herds within each grid cell on each day, and therefore a distribution of values of (R) for herds in each grid cell and on each day.
    We present the distribution of (R) for herds by considering the reproductive ratio that would be calculated if the random variable in Eq. (8), ({{b}_{l}}^{(h)}cdot {Z}_{lt}), took its pth percentile value every day, denoting this reproductive ratio, ({R}_{p}(x,t)). ({R}_{p}(x,t)) estimates the reproductive ratio that p% of herd reproductive ratios are less than in grid cell x on day t. Also, we numerically invert the threshold relationship to find the percentage value, ({varvec{P}}), such that ({R}_{p}(x,t)) satisfies the threshold quantity,

    $${varvec{P}}left(x,tright)={1-p in [mathrm{0,1}] | {R}_{p}(x,t) = 1 }.$$
    (9)

    ({varvec{P}}(x,t)) is therefore an estimate of the percentage of herds that could have a multiplying BTV outbreak in grid cell x if BTV was introduced on day t.
    To enable spatially extending risk predictions for BTV across Europe certain assumptions about how the midge biting model could be interpolated between the trap locations were necessary. The best regression model for midge catches found significantly different seasonality at the Italian catch locations compared to Sweden and The Netherlands. We assumed that day length was a determinant of midge seasonality and overwintering at different latitudes; for example, the first day of the year with a day length shorter than 9 h has been associated with the onset of overwintering in UK midges66. We found no midges on days with less than 8.5 h of daylight when trapping, although catching was not attempted outside of March-November so the number of samples was small. Therefore, at latitudes where no day is ever shorter than 8.5 h (lower latitudes than 46oN, which is more or less the border of Switzerland and Italy) we used Italian seasonality to predict midge biting. At latitudes that have at least one day shorter than 8 h (higher latitudes than 49oN) we used the Swedish and Dutch midge seasonality. Between these two latitudes, a linear interpolation between the predictions of the two seasonal models was applied.
    All map images were generated using MATLAB® contourf function. Both the livestock density and E-OBS datasets included grid cells that contain only water, these cells were coloured white.
    The MATLAB® code for generating the reproductive ratio estimates is available at a github repository (https://github.com/SamuelBrand1/BTV-European-Reproductive-Ratios). More