Ecology

1.EFSA. Effectiveness of in planta control measures for Xylella fastidiosa. EFSA J. 17(5). https://doi.org/10.2903/j.efsa.2019.5666 (2019).2.Hopkins, D. L. Xylella fastidiosa: Xylem-limited bacterial pathogen of plants. Annu. Rev. Phytopathol. 27(1), 271–290. https://doi.org/10.1146/annurev.py.27.090189.001415 (1989).Article 

Google Scholar 
3.Saponari, M., Boscia, D., Nigro, F. & Martelli, G. P. Identification of Dna sequences related to Xylella fastidiosa in oleander, almond and olive trees exhibiting leaf scorch symptoms in Apulia (southern Italy). J. Plant Pathol. 95(3), 668. https://doi.org/10.4454/JPP.V95I3.035 (2013).Article 

Google Scholar 
4.EPPO. Xylella fastidiosa in EPPO region. EPPO Bulletin. 49(2) (2019).5.Fierro, A., Liccardo, A. & Porcelli, F. A lattice model to manage the vector and the infection of the Xylella fastidiosa on olive trees. Sci. Rep. 9, 8723. https://doi.org/10.1038/s41598-019-44997-4 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
6.Saponari, M., Giampetruzzi, A., Loconsole, G., Boscia, D. & Saldarelli, P. Xylella fastidiosa in olive in apulia: Where we stand. Phytopathology 109(2), 175–186. https://doi.org/10.1094/PHYTO-08-18-0319-FI (2019).CAS 
Article 
PubMed 

Google Scholar 
7.Mannino, M. R. et al. Horizon scanning for plant health: Report on 2017–2020 activities. EFSA Support. Publ. https://doi.org/10.2903/sp.efsa.2021.EN-2010 (2021).Article 

Google Scholar 
8.EFSA. Scientific opinion on the risks to plant health posed by Xylella fastidiosa in the EU territory, with the identification and evaluation of risk reduction options. EFSA J. 13(1), 3989. https://doi.org/10.2903/j.efsa.2015.3989 (2015).9.Cornara, D. et al. An overview on the worldwide vectors of Xylella fastidiosa. Entomol. Gen. 39(3–4), 157–181. https://doi.org/10.1127/entomologia/2019/0811 (2019).Article 

Google Scholar 
10.Finke, D. L. Contrasting the consumptive and non-consumptive cascading effects of natural enemies on vector-borne pathogens. Entomol. Exp. Appl. 144, 45–55. https://doi.org/10.1111/j.1570-7458.2012.01258.x (2012).Article 

Google Scholar 
11.Martini, X., Hoffmann, M., Coy, M. R., Stelinski, L. L. & Pelz-Stelinski, K. S. Infection of an insect vector with a bacterial plant pathogen increases its propensity for dispersal. PLoS ONE 10(6), 1–16. https://doi.org/10.1371/journal.pone.0129373 (2015).CAS 
Article 

Google Scholar 
12.Almeida, R. P. P. et al. Addressing the new global threat of Xylella fastidiosa. Phytopathology 109(2), 172–174. https://doi.org/10.1094/PHYTO-12-18-0488-FI (2019).CAS 
Article 
PubMed 

Google Scholar 
13.Cornara, D., Bosco, D. & Fereres, A. Philaenus spumarius: When an old acquaintance becomes a new threat to European agriculture. J. Pest. Sci. 91(3), 957–972. https://doi.org/10.1007/s10340-018-0966-0 (2018).Article 

Google Scholar 
14.Halkka, O., Raatikainen, M., Vasarainen, A. & Heinonen, L. Ecology and ecological genetics of Philaenus spumarius (L.) (Homoptera). Ann. Zool. Fenn. 4, 1–18 (1967).
Google Scholar 
15.Lavigne, R. Biology of Philaenus leucophthalmus (L.) in Massachusetts. J. Econ. Entomol. 52(5), 904–907. https://doi.org/10.1093/jee/52.5.904 (1959).Article 

Google Scholar 
16.Ossiannilsson, F. The Auchenorrhyncha (Homoptera) of Fennoscandia and Denmark. Part 2: The families Cicadidae, Cercopidae, Membracidae, and Cicadellidae (excl. Deltocephalinae). Fauna Entomol. Scand. 7(2), 223–593 (1981).
Google Scholar 
17.Weaver, C. R. The seasonal behavior of meadow spittlebug and its relation to a control method. J. Econ. Entomol. 44(3), 350–353. https://doi.org/10.1093/jee/44.3.350 (1951).MathSciNet 
CAS 
Article 

Google Scholar 
18.Weaver, C. R. & King, D. R. Meadow spittlebug, Philaenus leucophthalmus (L.). Research Bulletin; Ohio Agricultural Experiment Station. (ed. Wooster, OH, USA, 1954).19.Drosopoulos, S. & Asche, M. Biosystematic studies on the spittlebug genus Philaenus with the description of a new species. Zool. J. Linn. Soc. 101(2), 169–177. https://doi.org/10.1111/j.1096-3642.1991.tb00891.x (2008).Article 

Google Scholar 
20.Grant, J. F., Lambdin, P. L. & Folium, R. A. Infestation levels and seasonal incidence of the meadow spittlebug (Homoptera: cercopidae) on musk thistle in Tennessee. J. Agric. Urban Entomol. 15, 83–91 (1998).
Google Scholar 
21.Halkka, O. Equilibrium populations of Philaenus spumarius L. Nature 193(4810), 93–94. https://doi.org/10.1038/193093a0 (1962).ADS 
Article 

Google Scholar 
22.Freeman, J. A. Studies in the distribution of insects by Aerial currents. J. Anim. Ecol. 14, 128 (1945).Article 

Google Scholar 
23.Reynolds, D. R., Chapman, J. W. & Stewart, A. J. A. Windborne migration of Auchenorrhyncha (Hemiptera) over Britain. Eur. J. Entomol. 114, 554–564. https://doi.org/10.14411/eje.2017.070 (2017).Article 

Google Scholar 
24.Gutierrez, A. P., Nix, H. A., Havenstein, D. E. & Moore, P. A. The ecology of Aphis Craccivora Koch and subterranean clover stunt virus in south-east Australia. III. A regional perspective of the phenology and migration of the Cowpea Aphid. J. Appl. Ecol. 11(1), 21–35. https://doi.org/10.2307/2402002 (1974).Article 

Google Scholar 
25.Pienkowski, R. L. & Medler, J. T. Synoptic weather conditions associated with long-range movement of the potato leafhopper, Empoasca fabae, into Wisconsin. Ann. Entomol. Soc. Am. 57(5), 588–591. https://doi.org/10.1093/aesa/57.5.588 (1964).Article 

Google Scholar 
26.Drake, V. A. Radar observations of moths migrating in a nocturnal low-level jet. Ecol. Entomol. 10(3), 259–265. https://doi.org/10.1111/j.1365-2311.1985.tb00722.x (1985).Article 

Google Scholar 
27.Wallin, J. R. & Loonan, D. V. Low-level jet winds, aphid vectors, local weather, and barley yellow dwarf virus outbreaks. Phytopathology 61(9), 1068. https://doi.org/10.1094/PHYTO-61-1068 (1971).Article 

Google Scholar 
28.Sedlacek, J. D. & Freytag, P. H. Aspects of the field biology of the Blackfaced Leafhopper (Homoptera: Cicadellidae) in corn and pastures in Kentucky. J. Econ. Entomol. 79(3), 605–613. https://doi.org/10.1093/jee/79.3.605 (1986).Article 

Google Scholar 
29.Zhu, M., Radcliffe, E. B., Ragsdale, D. W., MacRae, I. V. & Seeley, M. W. Low-level jet streams associated with spring aphid migration and current season spread of potato viruses in the U.S. northern Great Plains. Agric. For. Meteorol. 138(1–4), 192–202. https://doi.org/10.1016/j.agrformet.2006.05.001 (2006).ADS 
Article 

Google Scholar 
30.Bodino, N. et al. Dispersal of Philaenus spumarius (Hemiptera: Aphrophoridae), a vector of Xylella fastidiosa, in olive grove and meadow agroecosystems. Environ. Entomol. https://doi.org/10.1093/ee/nvaa140 (2020).Article 
PubMed Central 

Google Scholar 
31.Lago, C. et al. Dispersal of Neophilaenus campestris, a vector of Xylella fastidiosa, from olive groves to over-summering hosts. J. Appl. Entomol. https://doi.org/10.1111/jen.12888 (2021).Article 

Google Scholar 
32.Minter, M. et al. The tethered flight technique as a tool for studying life-history strategies associated with migration in insects. Ecol. Entomol. 43(4), 397–411. https://doi.org/10.1111/een.12521 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
33.Ávalos-Masó, J. A., Martí-Campoy, A. & Soto, T. A. Study of the flying ability of Rhynchophorus ferrugineus (Coleoptera: Dryophthoridae) adults using a computer-monitored flight mill. Bull. Entomol. Res. 104(4), 462–470. https://doi.org/10.1017/S0007485314000121 (2014).Article 

Google Scholar 
34.Yu, E. Y., Gassmann, A. J. & Sappington, T. W. Using flight mills to measure flight propensity and performance of western corn rootworm, diabrotica virgifera virgifera (Leconte). J. Vis. Exp. 152, e59196. https://doi.org/10.3791/59196 (2019).Article 

Google Scholar 
35.Riley, J. R., Downham, M. C. A. & Cooter, R. J. Comparison of the performance of Cicadulina leafhoppers on flight mills with that to be expected in free flight. Entomol. Exp. Appl. 83(3), 317–322. https://doi.org/10.1046/j.1570-7458.1997.00186.x (1997).Article 

Google Scholar 
36.Zhang, Y., Wang, L., Wu, K., Wyckhuys, K. A. G. & Heimpel, G. E. Flight performance of the Soybean Aphid, Aphis glycines (Hemiptera: Aphididae) under different temperature and humidity regimens. Environ. Entomol. 37(2), 301–306. https://doi.org/10.1603/0046-225X(2008)37[301:FPOTSA]2.0.CO;2 (2008).Article 
PubMed 

Google Scholar 
37.Cheng, Y., Luo, L., Jiang, X. & Sappington, T. Synchronized oviposition triggered by migratory flight intensifies larval outbreaks of beet. PLoS ONE https://doi.org/10.1371/journal.pone.0031562 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
38.Jones, C. M. et al. Genomewide transcriptional signatures of migratory flight activity in a globally invasive insect pest. Mol. Ecol. 24(19), 4901–4911. https://doi.org/10.1111/mec.13362 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
39.White, S. M., Bullock, J. M., Hooftman, D. A. P. & Chapman, D. S. Modelling the spread and control of Xylella fastidiosa in the early stages of invasion in Apulia, Italy. Biol. Invasions 19(6), 1825–1837. https://doi.org/10.1007/s10530-017-1393-5 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
40.Jones, V. P., Naranjo, S. E. & Smith, T. J. Insect ecology and behavior: Laboratory flight mill studies. Accessed 22 July 2021. (2010). http://entomology.tfrec.wsu.edu/VPJ_Lab/Flight-Mill41.Martí-Campoy, A. et al. Design of a computerised flight mill device to measure the flight potential of different insects. Sensors (Switzerland) 16(4), 485. https://doi.org/10.3390/s16040485 (2016).ADS 
Article 

Google Scholar 
42.Kees, A. M., Hefty, A. R., Venette, R. C., Seybold, S. J. & Aukema, B. H. Flight capacity of the walnut twig beetle (coleoptera: Scolytidae) on a laboratory flight mill. Environ. Entomol. 46(3), 633–641. https://doi.org/10.1093/ee/nvx055 (2017).Article 
PubMed 

Google Scholar 
43.Morente, M. et al. Distribution and relative abundance of insect vectors of Xylella fastidiosa in olive groves of the Iberian peninsula. Insects 9(4), 175. https://doi.org/10.3390/insects9040175 (2018).Article 
PubMed Central 

Google Scholar 
44.Morente, M., Cornara, D., Moreno, A. & Fereres, A. Continuous indoor rearing of Philaenus spumarius, the main European vector of Xylella fastidiosa. J. Appl. Entomol. 142(9), 901–904. https://doi.org/10.1111/jen.12553 (2018).Article 

Google Scholar 
45.Guthery, F. S., Burnham, K. P. & Anderson, D. R. Model Selection and multimodel inference: A practical information-theoretic approach. J. Wildl. Manag. 67, 655 (2003).Article 

Google Scholar 
46.Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A. & Smith, G. M. Mixed Effects Models and Extensions in Ecology with R. (ed. Springer Sci. Bus. Media, 2009).47.Strona, G., Carstens, C. J. & Beck, P. S. A. Network analysis reveals why Xylella fastidiosa will persist in Europe. Sci. Rep. 7(1), 1–8. https://doi.org/10.1038/s41598-017-00077-z (2017).CAS 
Article 

Google Scholar 
48.Whittaker, J. B. Density regulation in a population of Philaenus spumarius (L.) (Homoptera: Cercopidae). J. Anim. Ecol. 42(1), 163–172. https://doi.org/10.2307/3410 (1973).Article 

Google Scholar 
49.Wiman, N. G., Walton, V. M., Shearer, P. W., Rondon, S. I. & Lee, J. C. Factors affecting flight capacity of brown marmorated stink bug, Halyomorpha halys (Hemiptera: Pentatomidae). J. Pest Sci. 88(1), 37–47. https://doi.org/10.1007/s10340-014-0582-6 (2015).Article 

Google Scholar 
50.Strona, G. et al. Small world in the real world: Long distance dispersal governs epidemic dynamics in agricultural landscapes. Epidemics 30, 100384. https://doi.org/10.1016/j.epidem.2020.100384 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
51.Irwin, M. E. & Tresh, J. M. Long-range aerial dispersal of cereal aphids as virus vectors in North America. Philos. Trans. R. Soc. London. B Biol. Sci. 321(1207), 421–446. https://doi.org/10.1098/rstb.1988.0101 (1988).ADS 
Article 

Google Scholar 
52.Chapman, J. W., Reynolds, D. R. & Wilson, K. Long-range seasonal migration in insects: Mechanisms, evolutionary drivers and ecological consequences. Ecol. Lett. 18(3), 287–302. https://doi.org/10.1111/ele.12407 (2015).Article 
PubMed 

Google Scholar 
53.Fereres, A., Irwin, M. E. & Kampmeier, G. E. Aphid movement: Process and consecuences. in Aphids as crop pests. (ed.2 Emden, H. F. van, Harrington, R.). 196–224. https://doi.org/10.1079/9781780647098.0196 (CABI Publishing, 2017).54.Petrovskii, S., Mashanova, A. & Jansen, V. A. A. Variation in individual walking behavior creates the impression of a Levy flight. PNAS 108, 8704–8707. https://doi.org/10.1073/pnas.1015208108 (2011).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
55.Okano, K. Sublethal effects of a neonicotinoid insecticide on the sharpshooter vectors of Xylella fastidiosa. Doctoral dissertation (UC Berkeley, 2009).56.Robinet, C., David, G. & Jactel, H. Modeling the distances traveled by flying insects based on the combination of flight mill and mark-release-recapture experiments. Ecol. Modell. 402, 85–92. https://doi.org/10.1016/j.ecolmodel.2019.04.006 (2019).Article 

Google Scholar 
57.Taylor, R. A. J., Bauer, L. S., Poland, T. M. & Windell, K. N. Flight performance of agrilus planipennis (Coleoptera: Buprestidae) on a flight mill and in free flight. J. Insect Behav. 23(2), 128–148. https://doi.org/10.1007/s10905-010-9202-3 (2010).Article 

Google Scholar 
58.Srygley, R. B. & Lorch, P. D. Coping with uncertainty: Nutrient deficiencies motivate insect migration at a cost to immunity. Integr. Comp. Biol. 53, 1002–1013. https://doi.org/10.1093/icb/ict047 (2013).CAS 
Article 
PubMed 

Google Scholar 
59.Nilakhe, S. S. & Buainain, C. M. Observations on movement of spittlebug adults. Pesqui. Agropecuária Bras. Brasília 23, 123–134 (1988).
Google Scholar 
60.Neuman-Lee, L. A., Hopkins, G. R., Brodie, E. D. & French, S. S. Sublethal contaminant exposure alters behavior in a common insect: Important implications for trophic transfer. J. Environ. Sci. Heal. Part B Pestic. Food Contam. Wastes 48(6), 442–448. https://doi.org/10.1080/03601234.2013.761839 (2013).CAS 
Article 

Google Scholar 
61.Wilson, D. M. The central nervous control of flight in a locust. J. Exp. Biol. 38(2), 471–490 (1961).Article 

Google Scholar 
62.Yamanaka, T., Tatsuki, S. & Shimada, M. Flight characteristics and dispersal patterns of fall webworm (Lepidoptera: Arctiidae) males. Environ. Entomol. 30(6), 1150–1157. https://doi.org/10.1603/0046-225X-30.6.1150 (2001).Article 

Google Scholar 
63.Blackmer, J. L., Hagler, J. R., Simmons, G. S. & Henneberry, T. J. Dispersal of Homalodisca vitripennis (Homoptera: Cicacellidae) from a point release site in citrus. Environ. Entomol. 35(6), 1617–1625. https://doi.org/10.1093/ee/35.6.1617 (2006).Article 

Google Scholar 
64.Bodino, N. et al. Phenology, seasonal abundance and stage-structure of spittlebug (Hemiptera: Aphrophoridae) populations in olive groves in Italy. Sci. Rep. 9(1), 1–17. https://doi.org/10.1038/s41598-019-54279-8 (2019).CAS 
Article 

Google Scholar 
65.Minuz, R. L., Isidoro, N., Casavecchia, S., Burgio, G. & Riolo, P. Sex-dispersal differences of four phloem-feeding vectors and their relationship to wild-plant abundance in vineyard agroecosystems. J. Econ. Entomol. 106(6), 2296–2309. https://doi.org/10.1603/ec13244 (2013).CAS 
Article 
PubMed 

Google Scholar 
66.Waloff, N. Dispersal by flight of leafhoppers (Auchenorrhyncha: Homoptera). J. Appl. Ecol. 10, 705 (1973).Article 

Google Scholar 
67.Johnson, C. G. Physiological factors in insect migration by flight. Nature 198(4879), 423–427. https://doi.org/10.1038/198423a0 (1963).ADS 
Article 

Google Scholar 
68.Drake, V. A. & Gatehouse, A. G. Insect Migration. Tracking Resources through Space and Time. (ed. Cambridge University Press). 7(3) Cambridge UK. https://doi.org/10.1007/s10841-006-9039-4 (1995).69.Sappington, T. W. & Showers, W. B. Reproductive maturity, mating status, and long-duration flight behavior of agrotis ipsilon (Lepidoptera: Noctuidae) and the conceptual misuse of the oogenesis flight syndrome by entomologists. Environ. Entomol. 21(4), 677–688. https://doi.org/10.1093/ee/21.4.677 (1992).Article 

Google Scholar 
70.Zhao, X. C. et al. Does the onset of sexual maturation terminate the expression of migratory behaviour in moths? A study of the oriental armyworm, Mythimna separata. J Insect Physiol. 55(11), 1039–432009. https://doi.org/10.1016/j.jinsphys.2009.07.007 (2009).CAS 
Article 
PubMed 

Google Scholar 
71.Tigreros, N. & Davidowitz, G. Flight-fecundity tradeoffs in wing-monomorphic insects. Adv. Insect Phys. 56, 1–41. https://doi.org/10.1016/bs.aiip.2019.02.001 (2019).Article 

Google Scholar 
72.Drake, V. A. & Farrow, R. A. The influence of atmospheric structure and motions on insect migration. Ann. Rev. Entomol. 33(1), 183–210. https://doi.org/10.1146/annurev.en.33.010188.001151 (1988).Article 

Google Scholar 
73.Burt, P. J. A. & Pedgley, D. E. Nocturnal insect migration: Effects of local winds. Adv. Ecol. Res. 27, 61–92. https://doi.org/10.1016/S0065-2504(08)60006-9 (1997).Article 

Google Scholar 
74.Gordh, G. & McKirdy, S. The Handbook of Plant Biosecurity (Springer, 2014). https://doi.org/10.1007/978-94-007-7365-3Book 

Google Scholar  More

SubjectsIn total, twenty-two black-capped chickadees (nine males and 13 females) were tested between May and December 2019, and 16 black-capped chickadees (seven males and nine females) completed the experiment. One male and one female failed to learn Pretraining, and one female failed to learn Non-differential training (see descriptions below for training information); as a result, all three were removed from the experiment. In addition, one male and two females died of natural causes during the course of the study (see Ethics Declaration). For all birds, sex was determined by deoxyribonucleic acid analysis of blood samples37. All birds were captured in Edmonton (North Saskatchewan River Valley, 53.53°N, 113.53°W; Mill Creek Ravine, 53.52°N, 113.47°W), Alberta, Canada in January 2018 and January 2019 and were at least one year of age at capture, verified by examining outer tail rectrices38.Prior to the current experiment, all chickadees were individually housed in Jupiter Parakeet cages (30 × 40 × 40 cm; Rolf C. Hagen, Inc., Montreal, QB, Canada) in a single colony room. Therefore, birds did not have physical contact with each other, but did have visual and auditory contact. Birds had ad libitum access to food (Mazuri Small Bird Maintenance Diet; Mazuri, St. Louis, MO, USA), water with vitamins supplemented on alternating days (Monday, Wednesday, Friday; Prime Vitamin Supplement; Hagen, Inc.), a cup containing grit, and a cuttlebone. Additional nutritional supplements included three to five sunflower seeds daily, one superworm (Zophobas morio) three times a week, and a mixture of hard-boiled eggs and greens (spinach or parsley) twice a week. The colony rooms were maintained at approximately 20 °C and on a light:dark cycle that followed the natural light cycle for Edmonton, Alberta, Canada.One bird had previous experience with one operant experiment involving chick-a-dee calls but showed no difference in responding in comparison to the naïve birds. The remaining 15 birds had no previous experimental experience with black-capped chickadee-produced fee-bee songs or any experimental paradigm.Recordings of acoustic stimuliThe following acoustic stimuli were used in our previous published operant study which indicated that male and female chickadees can identify individual females via their song36. Stimuli included the songs of six female black-capped chickadees. All females were captured in Edmonton (North Saskatchewan River Valley, 53.53°N, 113.53°W; Mill Creek Ravine, 53.52°N, 113.47°W), Alberta, Canada in January 2010, 2011, 2012, and 2014, and all females were at least one year of age at capture, verified by examining outer tail rectrices38. Four females were recorded in Spring 2012 and two females were recorded in Fall 2014. Each recording session lasted approximately 1 h and all recordings took place after colony lights turned on at 08:00, specifically at 8:15. All females were recorded in silence, individually, within their respective colony room cages. Colony room cages were placed in sound-attenuating chambers for recording (1.7 m × 0.84 m × 0.58 m; Industrial Acoustics Company, Bronx, NY). An AKG C 1000S (AKG Acoustics, Vienna, Austria) microphone (positioned 0.1 m above and slightly behind the cage) was connected to a Marantz PMD670 (Marantz America, Mahwah, NJ) digital recorder (16 bit, 44,100 Hz sampling rate) and was used for all recordings. Audio recordings were analyzed and cut into individual files (songs) using SIGNAL 5.03.11 software (Engineering Design, Berkley, CA, USA).Acoustic stimuliFor the current study, a total of 150 vocalizations were used as stimuli, these vocalizations were comprised of 25 fee-bee songs produced by each of six recorded female chickadees. We ensured that all 150 were of high quality, meaning no audible interference, and all stimuli were bandpass filtered (lower bandpass 500 Hz, upper bandpass 14,000 Hz) using GoldWave version 6.31 (GoldWave, Inc., St. John’s, NL, Canada) in order to reduce any background noise outside of the song stimuli spectrum. For each song stimulus, 5 ms of silence was added to the leading and trailing portion of the vocalization and each stimulus was tapered to remove transients, in addition amplitude was equalized peak to peak using SIGNAL 5.03.11 software. When triggered, stimuli were presented at approximately 75 dB peak SPL as measured by a calibrated Brüel & Kjær Type 2239 (Brüel & Kjær Sound & Vibration Measurement A/S, Nærum, Denmark) sound pressure meter (A-weighting, slow response), a level that corresponds with the natural chickadee vocalizations amplitudes39,40,41. All dB measurements were made at the level of the request perch where birds trigger stimuli and where birds are required to remain for the length of the stimuli and all dB measurements refer to SPL.Noise stimuliAnthropogenic noise stimuli were originally created and used by Potvin and MacDougall-Shackleton42 and by Potvin, Curcio, Swaddle, and MacDougall-Shackleton43. The stimuli were recorded from an urban area in Melbourne, Victoria, Australia and other anthropogenic noise stimuli of various trains, cars, motorcycles, and lawnmowers downloaded from Soundbible.com were used. Within Victoria44 and Alberta45,46, urban traffic noise averages 60–80 dB SPL. The files used varied in length, with those recorded in Melbourne all being 10 min in length and those downloaded from Soundbible.com varying between 1–10 minutes42,43. In total 10 tracks were used with 30 total minutes of noise stimuli. Three anthropogenic noise conditions were used in the study, including Silence (no noise), Low noise (anthropogenic noise stimuli played at ~ 40 dB peak SPL), and High noise (anthropogenic noise stimuli played at ~ 75 dB peak SPL) replicating the variation of traffic noise experienced in urban areas42,43. For the Low and High noise conditions the 10 tracks repeated on a randomized loop during data collection (natural light of light/dark cycle) with, thus noise exemplars overlapped songs by chance, to further emulate urban areas. Noise stimuli had natural variations and modulations in frequency and amplitude over the course of the sound files. All dB measurements for noise stimuli included in this study refer to SPL. See Fig. 1 for female song and traffic noise stimuli spectrograms and power spectra.Figure 1(A) Spectrogram of a female fee-bee song in silence. (B) Power spectrum of female fee-bee song in silence . (C) Spectrogram of female fee-bee song in low noise. (D) Power spectrum of female fee-bee song (black) in low noise (grey). (E) Spectrogram of female fee-bee song in high noise. (F) Power spectrum of female fee-bee song (black) in high noise (grey).Full size imageApparatusFor the duration of the experiment, birds were housed individually in modified colony room cages (30 × 40 × 40 cm; described above) which were placed inside a ventilated, sound-attenuating operant chamber. See Fig. 2 for illustration of operant conditioning chamber. All chambers were lit with a full spectrum LED bulb (3 W, 250 lm E26, Not-Dim, 5000 K; Lohas LED, Chicago, IL, USA), and maintained the natural light:dark cycle for Edmonton, Alberta. Each cage within each operant chamber contained two perches and an additional perch fitted with an infrared sensor (i.e., the request perch). See Fig. 2C. Each cage also contained a water bottle, grit cup, and cuttlebone See Fig. 2G-2H. Birds had ad libitum access to water (with vitamins supplemented on alternating days; Monday, Wednesday, Friday), grit, and cuttlebone and were provided two superworms daily (a morning and afternoon worm). An opening (11 × 16 cm) located on the left side of the cage allowed the birds to access a motorized feeder, with a red LED light, and equipped with an infrared sensor47. See Fig. 2B,D–F. The purpose of the sensor was so that food was only available as a reward for correct responses to auditory stimuli during the operant discrimination task. We should note that performance of the discrimination task is required for access to food and thus maintains motivation. For operation and data collection, a personal computer connected to a single-board computer48 scheduled trials and recorded responses to stimuli. Stimuli were played from a personal computer hard drive through a Cambridge Integrated Amplifier (model A300 or Azur 640A; Cambridge Audio, London, England). Data is downloaded once a day in order to reduce stress on subjects as all equipment must be tested following download, requiring contact with subjects. Stimuli played in the chamber through a Fostex full-range speaker (model FE108 Σ or FE108E Σ; Fostex Corp., Japan; frequency response range 80–18,000 Hz) located beside the feeder. See Sturdy and Weisman49 for a detailed description of the apparatus. See Fig. 2 for an illustration of the operant conditioning chamber set-up.Figure 2Illustration of the operant conditioning chamber, including: (A) speaker, (B) automated feeder, (C) request perch fitted with infrared photo-beam assembly, (D) feeder cup, (E) electrical inputs, (F) red LED, (G) water bottle, (H) and cuttlebone. Also shown is the feeder opening, and additional perches. To simplify, the sketch the front and floor of the chamber, and the enclosure’s acoustic lining are not included.Full size imageProcedureOperant conditioningOur current operant conditioning go/no-go set-up is used to understand how birds perceive auditory stimuli. By training the birds to respond to particular stimuli and withhold responding to other stimuli we can compare responses to both types of stimuli. The go/no-go paradigm requires the birds to learn which stimuli require correct responses (go), providing reinforcement (food), and which stimuli require birds to withhold responding (no-go), resulting in the avoidance of punishment (lights out).The current study follows nine stages, after learning to use the operant conditioning set-up, birds then go through Non-differential training (stage 1) where they will be exposed to all stimuli that will be used in the experiment and to ensure that the birds respond to the stimuli equivalently. Then birds complete Discrimination training (stage 2) where birds on two categories of sounds. One category is rewarded, the other category is punished. Then the Discrimination-85 (stage 3) phase prepares birds for future trials where there is no reward nor punishment. After this point, birds will follow three series (Silence; Low; High) of Discrimination-85 with noise (stage 4, 6, 8) and a corresponding Probe with noise (stage 5, 7, 9), meaning that that each subject will repeated the two discrimination tasks three times with different noise conditions, with the order of noise conditions randomized among individuals. The detailed procedures for each stage are described in the following.Non-differential trainingThe purpose of Non-differential training is to engender a high level of responding on all trials, across all stimuli. Once a bird learned to use the request perch fitted with a sensor as well as learned to use the feeder to obtain food then Pretraining began. During Pretraining, birds were trained to respond to a 1 s tone (1,000 Hz) in order to receive access to food. Pretraining occurred over an approximately 15-day period in order to allow acclimatization to the chamber, feeder, and speaker. Following Pretraining was Non-differential training. During Non-differential training, birds received food for responding to all fee-bee song stimuli. All trials began when a bird landed on the request perch and remained on the perch for between 900–1100 ms, at which point a randomly-selected song stimulus played. Songs were presented in random order from trial to trial until all 150 stimuli had been triggered and played without replacement; once all 150 stimuli were played, a new random sequence initiated. In the event that the bird left the request perch during a stimulus presentation, the trial was deemed interrupted, and resulted in a 30 s lights out of the operant chamber. If the bird entered the feeder within 1 s after the stimulus (any stimulus) was played, it was given 1 s access to food, followed by a 30 s intertrial interval. If a bird remained on the request perch during the stimulus presentation and the 1 s following the completion of the stimulus, then the bird received a 60 s intertrial interval with the lights on. Birds continued on Non-differential training until they completed six 450-trial blocks at ≥ 60% responding on average to all stimuli, at least four 450-trial blocks at ≤ 3% difference in responding to future rewarded versus future unrewarded Discrimination stimuli, at least four 450-trial blocks at ≤ 3% difference in responding to future rewarded versus unrewarded Discrimination stimuli. Then following a day of free feed (during which birds had ad libitum access to a food cup) birds completed a second round of Non-differential training in which they completed at least one 450-trial block that met each of the above requirements. A 450-trial block consisted of the bird experiencing each of the 150 stimuli three times. For the current study the average time to complete Non-differential training ranged from 10 to 41 days (M = 21.43, SD = 9).Discrimination trainingDiscrimination training procedures included only 114 out of the 150 training stimuli that were previously presented in non-differential training, and responses to these stimuli were now differentially reinforced. Specifically, correct responses to half of the stimuli (“rewarded stimuli”, S+) were positively reinforced with 1 s access to food, and incorrect responses to the other half (“unrewarded stimuli”, S−) were instead punished with a 30-s intertrial interval of lights off within the operant chamber. In regard to criterion, Discrimination training continued until a bird completed six 342-trial blocks with a discrimination ratio between their respective S+ and S− of greater than 0.80 with the last two blocks being consecutive. For discrimination ratio calculations see Response Measures below.The current subjects were randomly assigned to either a True category discrimination group (n = 10) or Pseudo category discrimination group (n = 6). Furthermore, chickadees in the True category discrimination group were divided into two subgroups: (a) True 1 (n = 5; three females and two males) discriminated between 57 rewarded fee-bee songs produced by three individual female chickadees (S+) and 57 unrewarded fee-bee songs produced by another three individual female chickadees (S−); and (b) True 2 (n = 5; two females and three males) discriminated between the same songs with opposite rewards, properly, the 57 rewarded (S+) fee-bee songs were the S− from True 1 and the 57 unrewarded (S−) fee-bee songs were the S+ from True 1. For birds in the True category discrimination the average number of blocks completed per day for Discrimination training ranged from 2.4–4.4 blocks (3.3 ± 0.7 blocks).In similitude, the Pseudo category discrimination group was divided into two subgroups: (a) Pseudo 1 (n = 3; two females and one male) discriminated between 57 randomly-selected rewarded (S+) fee-bee songs and 57 randomly-selected unrewarded (S−) fee-bee songs; and (b) the second subgroup Pseudo 2 (n = 3; one female and two males) discriminated between the same songs with opposite rewards, meaning, the 57 rewarded (S+) fee-bee songs were the S− from Pseudo 1 and the 57 unrewarded (S−) fee-bee songs were the S+ from Pseudo 1 (S+) fee-bee songs and 57 randomly-selected unrewarded (S−) fee-bee songs. To explicate, the purpose of the two Pseudo groups was to include a control in which subjects are required to memorize each vocalization independent of the producer rather than be trained to categorize songs according to individual chickadees as the True groups have been. All birds remained in their respective groups (True 1 and 2; Pseudo 1 and 2) for the duration of the study. For birds in the Pseudo category discrimination the average number of blocks completed per day for Discrimination training ranged from 3.3–6.1 blocks (4.34 ± 1.2 blocks).Discrimination-85 phaseDiscrimination-85 was identical to the above Discrimination training except that rewarded songs were reinforced with a reduced probability, P = 0.85. Therefore, for 15% of trials when a rewarded stimulus was played and a bird correctly responded, no access to food was triggered. Instead, a 30 s lights on intertrial interval occurred. The change in reinforcement occurs in order to prepare birds for Probe trials in which novel song stimuli were neither rewarded with access to food nor unrewarded with a lights out, instead nothing occurs. Discrimination-85 continued until birds completed two consecutive 342-trial blocks with a discrimination ratio of at least 0.80.Discrimination-85 phase with noiseAll subjects followed three series (Silence; Low; High) of Discrimination-85 with noise and a corresponding Probe with noise and the order of noise stimuli was randomly-selected for each bird. Discrimination-85 with noise was identical to the Discrimination-85 phase except one of the three noise stimuli conditions (Silence; Low noise, 40 dB SPL; High noise, 75 dB SPL) was played over the song stimuli. The noise stimuli condition was randomly-selected for each bird. Each bird went through three series of Discrimination-85 with noise (Silence; Low; High) until reaching criteria: two consecutive 342-trial blocks with a discrimination ratio of at least 0.80. Here, we were interested in how the addition of noise would impact discrimination between rewarded and unrewarded female song stimuli.Probe phase with noiseFollowing each Discrimination-85 phase with noise was a corresponding Probe phase with noise. During Probe the reinforcement contingencies from Discrimination-85 were maintained. In addition to the 114 stimuli from Discrimination training, this stage included 12 novel fee-bee songs (i.e., Probe stimuli), two from each of the six individual females. For True groups, six of these novel songs were categorized as P + and the other six as P-, based on whether they were produced by the same birds as the S+ or the S− training stimuli. For Pseudo groups, the novel songs were not assigned to categories. For both groups, the 12 novel stimuli were neither rewarded (no food access) nor unrewarded (no lights out). The birds completed six 126-trial blocks in which the 114 familiar discrimination stimuli repeated once per block and the 12 probe sequences played once per block. In addition, one of the three noise stimuli conditions (Silence; Low noise, 40 dB SPL; High noise, 75 dB SPL) was played over the song stimuli, and each bird went through three series of Probe with noise (Silence; Low; High) which corresponded with the birds previous Discrimination-85 phase with noise condition. Thus, all birds completed all three Discrimination-85 phases with noise conditions followed by the corresponding Probe with noise conditions, and the order of noise stimuli condition was randomly-selected for each bird. In Probe phases we are interested if subjects can categorize novel stimuli to previously rewarded or unrewarded female birds.Response measuresFor each 342-block trial during training (Discrimination-85 with noise; Probe with noise), proportion response was calculated (R + /(N-I)): R + represents the number of trials in which the bird went to the feeder, N represents the total number of trials, and I represents the number of interrupted trials in which the bird left the perch before the entire stimulus played. For Discrimination training and the Discrimination-85 phase, a discrimination ratio was calculated by dividing the mean proportion response to all S+ stimuli by the mean proportion response to S+ stimuli plus the mean proportion response to S− stimuli. A discrimination ratio = 0.50 specifies equal response to rewarded (S+) and unrewarded (S−) stimuli, a discrimination ratio = 1.00 specifies a perfect discrimination between S+ and S− stimuli. We also collected data regarding the number of blocks and days per stage (Discrimination training; Discrimination-85 training with noise) in order to examine the latency of discrimination learning.Statistical analysesAll statistical analyses were conducted using SPSS (Version 20, Chicago, SPSS Inc.). In order to compare the number of trials needed to reach criterion and the discrimination ratios between True and Pseudo groups for Discrimination Training we conducted an analysis of variance (ANOVA). For Discrimination-85 with noise (Silence, Low noise, High noise), an ANOVA was conducted to compare the number of trials needed to reach criterion and the discrimination ratios between True and Pseudo groups. We also conducted post-hoc tests in order to reveal any sex differences between groups.And for Discrimination-85 with noise and Probes with noise repeated measures ANOVA was conducted to compare proportion response to training stimuli and probe stimuli between True groups and Pseudo groups. Lastly, we conducted post-hoc tests in order to reveal any differences in the number of trials to reach criterion during Discrimination training and to Discrimination-85 with noise.Ethics declarationThroughout the experiment, birds remained in the testing apparatus to minimize the transport and handling of each bird. One male and two female subjects died from natural causes during operant training. Following the experiment, healthy birds were returned to the colony room for use in future experiments.All procedures were conducted in accordance with the Canadian Council on Animal Care (CCAC) Guidelines and Policies with approval from the Animal Care and Use Committee for Biosciences for the University of Alberta (AUP 1937), which is consistent with the Animal Care Committee Guidelines for the Use of Animals in Research. Birds were captured and research was conducted under an Environment Canada Canadian Wildlife Service Scientific permit (#13-AB-SC004), Alberta Fish and Wildlife Capture and Research permits (#56,066 and #56,065), and the City of Edmonton Parks permit. All methods are reported in accordance with ARRIVE guidelines. More

Resistance of H. zea populations to Cry1AcTo measure variation in resistance of H. zea populations across field locations during 2017 and 2018, larval offspring of insects collected from non-Bt maize were subjected to a diet-overlay bioassay containing a diagnostic concentration of Cry1Ac (29 µg/cm2) corresponding to the mean LC95 of four Cry1Ac susceptible H. zea populations. Overall, larval survivorship varied significantly among years (Fig. 1).Figure 1Survival of H. zea larvae collected in 2017 and 2018 following exposure to 29 µg/cm2 Cry1Ac in diet-overlay assay. Dashed line is mean survival of a known susceptible field population. ***Years significantly different F = 25.38; df = 1, 58; P  More

1.Lengeler, C. Insecticide-treated bed nets and curtains for preventing malaria. Cochrane Database Syst Rev. https://doi.org/10.1002/14651858.CD000363.pub2 (2004).Article 
PubMed 
PubMed Central 

Google Scholar 
2.Pluess, B. et al. Indoor residual spraying for preventing malaria (Review). Cochrane Rev. https://doi.org/10.1002/14651858.CD006657.pub2 (2010).Article 

Google Scholar 
3.WHO. Guidelines for malaria vector control. (World Health Organization, 2019). https://www.who.int/malaria/publications/atoz/9789241550499/en/. Accessed 08 Sept 2019.4.Bhatt, S. et al. The effect of malaria control on Plasmodium falciparum in Africa between 2000 and 2015. Nature 526(7572), 207–211 (2015).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Killeen, G. F. et al. Made-to-measure malaria vector control strategies: Rational design based on insecticide properties and coverage of blood resources for mosquitoes. Malar. J. 13(1), 146 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
6.Waite, J. L. et al. Increasing the potential for malaria elimination by targeting zoophilic vectors. Sci. Rep. 7, 40551 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
7.Cooke, M. K. et al. A bite before bed’: Exposure to malaria vectors outside the times of net use in the highlands of western Kenya. Malar. J. 14(1), 259 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
8.Milali, M. P., Sikulu-Lord, M. T. & Govella, N. J. Bites before and after bedtime can carry a high risk of human malaria infection. Malar. J. 16(1), 91 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
9.Pates, H. et al. Differential behaviour of Anopheles gambiae sensu stricto (Diptera: Culicidae) to human and cow odours in the laboratory. Bull. Entomol. Res. 91(4), 289–296 (2001).CAS 
PubMed 
Article 

Google Scholar 
10.Costantini, C. et al. Odor-mediated host preferences of West African mosquitoes, with particular reference to malaria vectors. Am. J. Trop. Med. Hyg. 58(1), 56–63 (1998).CAS 
PubMed 
Article 

Google Scholar 
11.Gillies, M. T. & DeMeillon, B. The Anophelinae of Africa South of the Sahara (Ethiopian Zoogeographical Region) (South African Institute for Medical Research, 1968).
Google Scholar 
12.Hamon, J. Les moustiques anthropophiles de la région de Bobo-Dioulasso (République de Haute-Volta): Cycles d’agressivité et variations saisonnières. Ann. Soc. Entomol. (1963).13.Nyarango, P. M. et al. A steep decline of malaria morbidity and mortality trends in Eritrea between 2000 and 2004: the effect of combination of control methods. Malar. J. 5(1), 33 (2006).PubMed 
PubMed Central 
Article 

Google Scholar 
14.WHO, World malaria report 2019, in Licence: CC BY-NC-SA 3.0 IGO. https://www.who.int/news-room/feature-stories/detail/world-malaria-report-2019. Accessed 17 Dec 2019.15.Hancock, P. A. et al. Mapping trends in insecticide resistance phenotypes in African malaria vectors. PLoS Biol. 18(6), e3000633 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Martinez-Torres, D. et al. Molecular characterization of pyrethroid knockdown resistance (kdr) in the major malaria vector Anopheles gambiae ss. Insect Mol. Biol. 7(2), 179–184 (1998).CAS 
PubMed 
Article 

Google Scholar 
17.Balabanidou, V. et al. Cytochrome P450 associated with insecticide resistance catalyzes cuticular hydrocarbon production in Anopheles gambiae. Proc. Natl. Acad. Sci. USA 113(33), 9268–9273 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Edi, C. V. et al. CYP6 P450 enzymes and ACE-1 duplication produce extreme and multiple insecticide resistance in the malaria mosquito Anopheles gambiae. PloS Genet 10(3), e1004236 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
19.Kreppel, K. et al. Emergence of behavioural avoidance strategies of malaria vectors in areas of high LLIN coverage in Tanzania. Sci. Rep. 10(1), 1–11 (2020).Article 
CAS 

Google Scholar 
20.Moiroux, N. et al. Human exposure to early morning Anopheles funestus biting behavior and personal protection provided by long-lasting insecticidal nets. PLoS ONE 9, e104967 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
21.Reddy, M. R. et al. Outdoor host seeking behaviour of Anopheles gambiae mosquitoes following initiation of malaria vector control on Bioko Island, Equatorial Guinea. Malar. J. 10(1), 184–184 (2011).PubMed 
PubMed Central 
Article 

Google Scholar 
22.Rozendaal, J. et al. Behavioral responses of Anopheles darlingi in Suriname to DDT residues on house walls. J. Am. Mosq. Control Assoc. 5(3), 339–350 (1989).CAS 
PubMed 

Google Scholar 
23.Mwangangi, J. M. et al. Shifts in malaria vector species composition and transmission dynamics along the Kenyan coast over the past 20 years. Malar. J. 12(1), 1 (2013).Article 

Google Scholar 
24.Sougoufara, S. et al. Shift in species composition in the Anopheles gambiae complex after implementation of long-lasting insecticidal nets in Dielmo, Senegal. Med. Vet. Entomol. 30(3), 365–368 (2016).CAS 
PubMed 
Article 

Google Scholar 
25.Moiroux, N. et al. Changes in Anopheles funestus biting behaviour following universal coverage of long-lasting insecticidal nets in Benin. J. Infect. Dis. 206, 1622–1629 (2012).CAS 
PubMed 
Article 

Google Scholar 
26.Matowo, N. S. et al. Using a new odour-baited device to explore options for luring and killing outdoor-biting malaria vectors: A report on design and field evaluation of the Mosquito Landing Box. Parasit. Vectors 6(1), 1–16 (2013).Article 

Google Scholar 
27.Russell, T. L. et al. Increased proportions of outdoor feeding among residual malaria vector populations following increased use of insecticide-treated nets in rural Tanzania. Malar. J. 10, 1 (2011).Article 

Google Scholar 
28.Sherrard-Smith, E. et al. Mosquito feeding behavior and how it influences residual malaria transmission across Africa. Proc. Natl. Acad. Sci. 116(30), 15086–15095 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
29.Dambach, P. et al. Reduction of malaria vector mosquitoes in a large-scale intervention trial in rural Burkina Faso using Bti based larval source management. Malar. J. 18(1), 311 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
30.Bayili, K. et al. Evaluation of efficacy of Interceptor (R) G2, a long-lasting insecticide net coated with a mixture of chlorfenapyr and alpha-cypermethrin, against pyrethroid resistant Anopheles gambiae s.l. in Burkina Faso. Malar. J. 16, 9 (2017).Article 
CAS 

Google Scholar 
31.Oumbouke, W. A. et al. Evaluation of an alpha-cypermethrin + PBO mixture long-lasting insecticidal net VEERALIN® LN against pyrethroid resistant Anopheles gambiae s.s.: an experimental hut trial in M’bé, central Côte d’Ivoire. Parasit. Vectors 12(1), 544 (2019).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
32.Yaro, J. B. et al. A cohort study to identify risk factors for Plasmodium falciparum infection in Burkinabe children: Implications for other high burden high impact countries. Malar. J. 19(1), 371 (2020).MathSciNet 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.Dabire, K. et al. Anopheles funestus (Diptera: Culicidae) in a humid savannah area of western Burkina Faso: Bionomics, insecticide resistance status, and role in malaria transmission. J. Med. Entomol. 44(6), 990–997 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
34.Badolo, A. et al. Three years of insecticide resistance monitoring in Anopheles gambiae in Burkina Faso: Resistance on the rise?. Malar. J. 11(1), 1 (2012).Article 

Google Scholar 
35.Hamon, J. et al. Présence dans le Sud-Ouest de la Haute-Volta d’une Population d’Anopheles gambiae” A” résistante au DDT (Organisation mondiale de la Santé, 1968).
Google Scholar 
36.Hamon, J. et al. Présence dans le Sud-Ouest de la Haute-Volta de populations d’Anopheles funestus Giles résistantes à la dieldrine. Med. Trop. 28(2), 222–226 (1968).
Google Scholar 
37.Toe, H. Characterisation of Insecticide Resistance in Anopheles gambiae from Burkina Faso and Its Impact on Current Malaria Control Strategies (University of Liverpool, 2015).
Google Scholar 
38.Namountougou, M. et al. Insecticide resistance mechanisms in Anopheles gambiae complex populations from Burkina Faso, West Africa. Acta Trop. 197, 105054 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Toe, K. H. et al. Increased pyrethroid resistance in malaria vectors and decreased bed net effectiveness, Burkina Faso. Emerg. Infect. Dis. 20, 10 (2014).Article 
CAS 

Google Scholar 
40.Beier, J. C., Killeen, G. & Githure, J. I. Entomologic inoculation rates and Plasmodium falciparum malaria prevalence in Africa. Am. J. Trop. Med. Hyg. 61, 109–113 (1999).CAS 
PubMed 
Article 

Google Scholar 
41.Smith, D. L. et al. The entomological inoculation rate and Plasmodium falciparum infection in African children. Nature 438, 492–495 (2005).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
42.Killeen, G. F. et al. Quantifying behavioural interactions between humans and mosquitoes: Evaluating the protective efficacy of insecticidal nets against malaria transmission in rural Tanzania. BMC Infect. Dis. 6, 161 (2006).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
43.Huho, B. et al. Consistently high estimates for the proportion of human exposure to malaria vector populations occurring indoors in rural Africa. Int. J. Epidemiol. 42(1), 235 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
44.Robert, V. & Carnevale, P. Influence of deltamethrin treatment of bed nets on malaria transmission in the Kou valley, Burkina Faso. Bull. World Health Organ. 69(6), 735 (1991).CAS 
PubMed 
PubMed Central 

Google Scholar 
45.Dabire, K. et al. Year to year and seasonal variations in vector bionomics and malaria transmission in a humid savannah village in west Burkina Faso. J. Vector Ecol. 33(1), 70–75 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Epopa, P. S. et al. Seasonal malaria vector and transmission dynamics in western Burkina Faso. Malar. J. 18(1), 113 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
47.Dambach, P. et al. Nightly biting cycles of anopheles species in rural northwestern Burkina Faso. J. Med. Entomol. 55(4), 1027–1034 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
48.Guglielmo, F. et al. Quantifying variation in exposure risk to mosquito bites at the individual level in Burkina Faso. Malar. J. (2020)49.Edi, C. A. et al. Long-term trends in Anopheles gambiae insecticide resistance in Côte d’Ivoire. Parasit. Vectors 7(1), 500 (2014).PubMed 
PubMed Central 

Google Scholar 
50.Fuseini, G. et al. Evaluation of the residual effectiveness of Fludora fusion WP-SB, a combination of clothianidin and deltamethrin, for the control of pyrethroid-resistant malaria vectors on Bioko Island, Equatorial Guinea. Acta Trop. 196, 42–47 (2019).CAS 
PubMed 
Article 

Google Scholar 
51.Rongsriyam, Y. & Busvine, J. Cross-resistance in DDT-resistant strains of various mosquitoes (Diptera, Culicidae). Bull. Entomol. Res. 65(3), 459–471 (1975).Article 

Google Scholar 
52.Chandre, F. et al. Status of pyrethroid resistance in Anopheles gambiae sensu lato. Bull. World Health Organ. 77(3), 230–234 (1999).CAS 
PubMed 
PubMed Central 

Google Scholar 
53.Bagi, J. et al. When a discriminating dose assay is not enough: measuring the intensity of insecticide resistance in malaria vectors. Malar. J. 14(1), 210 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
54.Robert, V., et al., Etude des taux de parturité et d’infection du complexe Anopheles gambiae dans la rizière de la vallée du Kou, Burkina Faso. Le paludisme en Afrique de l’Ouest, 1991: p. 17.55.Pombi, M. et al. The Sticky Resting Box, a new tool for studying resting behaviour of Afrotropical malaria vectors. Parasit. Vectors 7(1), 247 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
56.Gnémé, A. et al. Anopheline occurrence and the risk of urban malaria in the city of Ouagadougou, Burkina Faso. (2019).57.Akogbéto, M. C. et al. Blood feeding behaviour comparison and contribution of Anopheles coluzzii and Anopheles gambiae, two sibling species living in sympatry, to malaria transmission in Alibori and Donga region, northern Benin, West Africa. Malar. J. 17(1), 307 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
58.Wamae, P. M. et al. Early biting of the Anopheles gambiae s.s. and its challenges to vector control using insecticide treated nets in western Kenya highlands. Acta Trop. 150, 136–142 (2015).CAS 
PubMed 
Article 

Google Scholar 
59.Gillies, M. Age-groups and the biting cycle in Anopheles gambiae: A preliminary investigation. Bull. Entomol. Res. 48(3), 553–559 (1957).Article 

Google Scholar 
60.Bayoh, M. et al. Persistently high estimates of late night, indoor exposure to malaria vectors despite high coverage of insecticide treated nets. Parasit. Vectors 7, 380 (2014).PubMed 
Article 

Google Scholar 
61.Sougoufara, S. et al. Biting by Anopheles funestus in broad daylight after use of long-lasting insecticidal nets: A new challenge to malaria elimination. Malar. J. 13, 125 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
62.Meyers, J. I. et al. Increasing outdoor host-seeking in Anopheles gambiae over 6 years of vector control on Bioko Island. Malar. J. 15(1), 1 (2016).MathSciNet 
Article 

Google Scholar 
63.Bayoh, M. N. et al. Anopheles gambiae: Historical population decline associated with regional distribution of insecticide-treated bed nets in western Nyanza Province, Kenya. Malar. J. 9(1), 62 (2010).PubMed 
PubMed Central 
Article 

Google Scholar 
64.Chinula, D. et al. Proportional decline of Anopheles quadriannulatus and increased contribution of An. arabiensis to the An. gambiae complex following introduction of indoor residual spraying with pirimiphos-methyl: An observational, retrospective secondary analysis of pre-existing data from south-east Zambia. Parasit. Vectors 11(1), 544 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
65.Main, B. J. et al. The genetic basis of host choice and resting behavior in the major African malaria vector Anopheles arabiensis. BioRxiv 2016, 044701 (2016).
Google Scholar 
66.Coluzzi, M. et al. Behavioral divergences between mosquitos with different inversion karyotypes in polymorphic populations of Anopheles gambiae complex. Nature 266, 832–833 (1977).ADS 
CAS 
PubMed 
Article 

Google Scholar 
67.Elanga-Ndille, E. et al. The G119S Acetylcholinesterase (Ace-1) target site mutation confers carbamate resistance in the major malaria vector Anopheles gambiae from cameroon: A challenge for the coming IRS implementation. Genes 10(10), 790 (2019).CAS 
PubMed Central 
Article 
PubMed 

Google Scholar 
68.Toé, K. H. et al. The recent escalation in strength of pyrethroid resistance in Anopheles coluzzi in West Africa is linked to increased expression of multiple gene families. BMC Genomics 16(1), 146 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
69.Finda, M. F. et al. Linking human behaviours and malaria vector biting risk in south-eastern Tanzania. PLoS ONE 14(6), e0217414 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
70.INSD. Annuaire sttistique 2018 Burkina Faso: Institut national de la statistique et de la démographie (2019).71.QGIS.org. QGIS Geographic Information System. QGIS Association. http://www.qgis.org. (2021).72.OpenStreetMap. OpenStreetMap contributors. https://www.openstreetmap.org. (2021).73.Benedict, M. Q. Care and maintenance of anopheline mosquito colonies. In The Molecular Biology of Insect Disease Vectors 3–12 (Springer, 1997).Chapter 

Google Scholar 
74.WHO. Test procedures for insecticide resistance monitoring in malaria vector mosquitoes: 2nd edition. (World Health Organization, 2016). http://www.who.int/malaria.75.WHO. Supplies for Monitoring Insecticide Resistance in Disease Vectors. (World Health Organization: Parasitic Diseases and Vector Control (PVC)/Communicable Disease Control, Prevention and Eradication (CPE), 2001).76.M.W. Service. Mosquito Ecology Field Sampling Methods (Elsevier Science Publishers, 1993).Book 

Google Scholar 
77.Kreppel, K. S. et al. Comparative evaluation of the Sticky-Resting-Box-Trap, the standardised resting-bucket-trap and indoor aspiration for sampling malaria vectors. Parasit. Vectors 8(1), 1–5 (2015).Article 

Google Scholar 
78.Gillies, M. & Coetzee, M. A Supplement to the Anophelinae of Africa south of the Sahara (Afrotropical region). (1987).79.Fanello, C., Santolamazza, F. D. & DellaTorre, A. Simultaneous identification of species and molecular forms of the Anopheles gambiae complex by PCR-RFLP. Medi. Vet. Entomol. 16(4), 461–464 (2002).CAS 
Article 

Google Scholar 
80.Sanou, A. et al. Evaluation of mosquito electrocuting traps as a safe alternative to the human landing catch for measuring human exposure to malaria vectors in Burkina Faso. Malar. J. 18(1), 386 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
81.Wirtz, R. et al. ELISA method for detecting Plasmodium falciparum circumsporozoite antibody. Bull. World Health Organ. 67(5), 535 (1989).CAS 
PubMed 
PubMed Central 

Google Scholar 
82.WHO. Manual on Practical Entomology in Malaria. Part II. Methods and Techniques. 84–85 (WHO Division of Malaria and other Parasitic Diseases, 1975).83.Garrett-Jones, C. The human blood index of malaria vectors in relation to epidemiological assessment. Bull. World Health Organ. 30(2), 241 (1964).CAS 
PubMed 
PubMed Central 

Google Scholar 
84.Beier, J. C. et al. Bloodmeal identification by direct enzyme-linked immunosorbent assay (ELISA), tested on Anopheles (Diptera: Culicidae) in Kenya. J. Med. Entomol. 25(1), 9–16 (1988).ADS 
CAS 
PubMed 
Article 

Google Scholar 
85.Kiszewski, A. et al. A global index representing the stability of malaria transmission. Am. J. Trop. Med. Hyg. 70(5), 486–498 (2004).PubMed 
Article 

Google Scholar 
86.Team, R.C. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2018). https://www.R-project.org/.87.Wood, S. & Wood, M.S. The mgcv package. https://www.R-project.org/ (2007).88.Bates, D. et al. Fitting linear mixed-effects models using lme4. http://arxiv.org/abs/1406.5823 (2014). More

Isolation of marine sediment-derived actinomycetes from west-central PhilippinesThe geographical sites identified in this study were evaluated to explore the actinomycete distribution in west-central Philippines (Fig. 1a). A total of 16 sediment cores were recovered from the 11 sampling sites and were processed in the laboratory using cultured-dependent actinomycete isolation (Supplementary Table S2). The seawater’s physicochemical conditions in all sampling sites were uniform with pH 7.0 and salinity ranging from 3.1 to 3.2. The characteristics of marine sediments and the distance of the actual collection sites identified using the given criteria varied per sampling location. The sediment characteristics vary from coarse to fine sand with mixture of broken corals and pebbles (Supplementary Table S3). Culture-dependent isolation revealed that actinomycete-like colonies and spores were observed in the minimal marine media after 30 to 60 days of incubation. Actinomycete isolates were repeatedly sub-cultured in enriched marine medium 1 (MM1) to obtain pure cultures as shown in Fig. 1b. Actinomycete growth observed in enriched media was white to gray aerial spores with brown to yellow mycelia or without diffusible pigmentations. Notably, there were strains that produced colonies with no diffusible pigmentations.Figure 1Distribution, abundance, and antibacterial activity of marine sediment-derived actinomycetes in the Philippines. (a) The overall map showing the 11 sampling sites situated within the west-central region in the Philippines. The enlarged map showed the details of the sampling sites and their corresponding actinomycetes abundance. Sampling sites are number-coded as shown inside the circle. The color gradient indicates strain abundance. (b) Actinomycete colonies were sub-cultured repeatedly to obtain pure culture of isolates. (c) A total of 92 out of 2212 actinomycetes strains have confirmed antibacterial activities as verified using microbroth susceptibility assay. The map with Streptomyces abundance plot was generated using ggplot2 package in Rstudio ver. 1.2.5042 (https://www.rstudio.com/).Full size imageIn this work, a total of 2212 pure actinomycete strains were isolated from marine sediments collected in 11 geographically distant sampling sites across the west-central Philippines (Fig. 1a). Actinomycete strains were highly abundant in Negros Occidental with 580 isolates (26%), followed by Southern Antique with 348 isolates (16%) and Southeastern Iloilo with 228 (10%). We recovered least actinomycete strains in Occidental Mindoro and South Central Visayas with only 94 and 81 strains, respectively.Antibacterial activity profile of actinomycete strainsWe assessed the antibacterial activities of actinomycete strains against a multidrug-resistant Gram-positive bacterium (S. aureus ATCC BAA-44) and three Gram-negative bacteria (E. coli ATCC 25922, P. aeruginosa ATCC 27853, and E. aerogenes ATCC 13048) using resazurin agar overlay assay and microbroth susceptibility assay as initial and confirmatory screenings, respectively. A total of 218 (9.9%) out of the 2212 actinomycete isolates have antibacterial activities in the initial screening as indicated by positive results or retained blue resazurin color in wells containing actinomycete overlaid with the test pathogens (Supplementary Fig. S2). The 218 active isolates were fermented to produce biomass for secondary metabolite extraction and to confirm their antibacterial activities by microbroth susceptibility testing. Figure 1c showed the confirmed 92 (4.1%) antibiotic-producing actinomycete strains. The majority of the strains (71 isolates) exhibited activity against Gram-positive S. aureus ATCC BAA-44. Thirty-nine (39) strains (42%) were active against E. coli ATCC 25922. Six strains were active against P. aeruginosa ATCC 27853, while all strains tested were inactive against E. aerogenes ATCC 13048 as indicated with no or less than 50% growth inhibition. Twenty-three (23) active strains targeted 2–3 test pathogens, while 69 active strains were only active against one test pathogen (Supplementary Table S4).Phylogenetic diversity of multiple antibiotic-producing strainsThe 92 active actinomycete strains were further identified and confirmed as Streptomyces species based on genomic analysis of their 16S rRNA and rpoB gene sequences. Comparison of 16S rRNA gene sequences (ranging from 1150 to 1500 nucleotides) and rpoB (700–995 nucleotides) gene sequences with their similar matches in the GenBank verified that the 92 active strains were closely related (97 to 100%) with 19 species under the genus Streptomyces. The nearly complete 16S rRNA and rpoB gene sequences were analyzed in a phylogenetic tree using maximum likelihood algorithms. The 16S rRNA and rpoB gene sequences of active strains reported in the present study were deposited in the GenBank nucleotide database (Supplementary Table S5).Phylogenetic analysis revealed multiple strains with identical 16S rRNA gene sequences which clustered together into 13 major clusters (shown by the colored nodes in the tree) with high bootstrap values ( > 90%) in the phylogenetic tree (Supplementary Fig. S3). Thirty-three strains (36%) were highly similar to S. parvulus presented in red circle. Followed by 12 strains with high similarity to S. enissocaesilis (light blue circle), 11 S. rochei strains (dark pink), six S. mutabilis strains (dark blue), five S. diastaticus strains (light green), four S. kunmingensis strains (green) and three S. geysiriensis strains (light orange). A phylogenetic analysis of rpoB gene sequence was conducted to provide a better resolution of the evolutionary relationship among strains within and between species supporting the taxonomic identity of the phylogenetically identical strains.Similarly, the rpoB gene sequences phylogenetic tree (Fig. 2) showed 13 major clusters that were highly supported with bootstrap replicates  > 90%, except for monophyletic clusters III (Streptomyces sp. strain DSD176) and X (Streptomyces sp. strain DSD1006) with low bootstrap replicates ( 1% abundance) can be recovered in both methods. We also noted that some species were exclusively recovered using a specific method. Eight species were exclusively recovered in DSM compared to 6 species in HSM. All active strains of S. kunmingensis, S. mutabilis, S. sedi and S. olivaceus were only recovered by HSM. Contrary, the active strains of S. carpaticus and S. harbinensis were only isolated using DSM.Carbon source composition of marine minimal media, along with effective isolation techniques, was crucial for the isolation of antibiotic-producing Streptomyces. Among the five minimal marine media used, three carbon sources yielded high isolation rate: glucose, mannitol, and trehalose yielded nine species with 32 strains (35%), nine species with 18 strains (20%), and ten species with 18 strains (20%) respectively. However, only eight species (12 strains, 20%) and two species (12 strains, 20%) were isolated in raffinose, and starch-based media, respectively. In the contrary, high diversity indices were observed in trehalose (2.197), mannitol (2.0), and raffinose (1.979) (Fig. 5a). As expected, the starch-containing media had the lowest isolation rate and diversity as only two species (S. enissocaesilis and S. parvulus) were able to utilize a more complex carbon source.Figure 5Diversity of antibiotic-producing Streptomyces using five different carbon sources. (a) From the five carbon sources in the minimal marine media utilized by Streptomyces strains in this study, mannitol yielded the highest number of active strains, while high diversity was recorded in active strains that utilized glucose (n = 92). (b) Venn diagram of five carbon sources showed that two Streptomyces species can be isolated using all five carbon sources.Full size imageInterestingly, co-isolation of species in different carbon sources was shown in the Venn diagram (Fig. 5b). Eight species can be recovered from at least two different carbon sources, whereas 11 species were exclusively isolated from a specific carbon source. Bioactive S. enissocaesilis and S. parvulus strains were recovered from all of the carbon sources utilized in this study. Active S. rochei were isolated in four media but not in starch-based media. In contrast, more exclusive species were isolated in trehalose with four species (S. harbinensis strains, Streptomyces sp. strain DSD3025, S. pseudogriseolus, and S. xiamenensis). Followed by glucose with three species (Streptomyces sp. strain DSD742, S. carpaticus and S. sedi), and two species each on mannitol- and raffinose-based media (Fig. 5b). The results indicated that diverse Streptomyces species preferred simple sugars-containing one or two sugar molecule as nutrient source compared to complex sugars.
Streptomyces abundance and diversity in geographical sampling locationsBioactive Streptomyces species were widely distributed across the different sampling locations in west-central Philippines (Fig. 6a). Although Southern Antique, Negros Occidental and Negros Oriental have highest number of active strains isolated, we found that Southern Antique, Southern Iloilo, and Western Antique were the most diverse sampling sites (Fig. 6b). We have isolated the greatest number of antibiotic-producing Streptomyces species which were evenly distributed in Southern Antique. This indicates that Southern Antique is stable with many potential niches that can support highly diverse Streptomyces species.Figure 6(a) The abundance profile of 19 antibiotic-producing Streptomyces species in different geographical location across the Philippines revealed that S. parvulus was the most abundant species. (b) Diversity, richness and evenness of antibiotic-producing Streptomyces species in different sites were calculated using vegan package in RStudio ver1.2.5042 (https://www.rstudio.com/). (c) Recovery profile and abundance of antibiotic-producing strain per depth layer, indicated by colored circles and its size, revealed that strains in the bottom sediments were the most diverse as compared to other depth layer and S. parvulus was the most abundant species recovered.Full size imageStreptomyces parvulus emerged as the dominant antibiotic-producing species in this study. Out of the 33 active S. parvulus strains, the majority were isolated in Negros Occidental (11 strains) and Negros Oriental (ten strains). Although Negros Occidental and Negros Oriental have high abundance, its microbial community is highly dominated by one species, S. parvulus, supported by the low species richness, evenness, and diversity (Fig. 6b). This finding implies that these sites may have few potential niches that only a few species dominate. Bioactive S. enissocaesilis strains were recovered in four sampling locations only; specifically, Southern Antique (five strains), Western Antique (three strains), Southeastern Iloilo (three strains), and Northwestern Antique (one strain). Active S. rochei were isolated and evenly distributed in seven sampling locations, but were not present in Western Antique, Northwestern Antique, South Central Visayas, and Tubbataha Reefs. Notably, we observed that no bioactive S. parvulus, S. enissocaesilis, and S. rochei were isolated in Tubbataha Reefs, but antibiotic-producing S. cacaoi, S. psuedogriseolus and S. mutabilis strains were isolated only in Tubbataha Reefs marine sediments. Meanwhile, site-specific species such as S. sedi were recovered only in Occidental Mindoro. The isolation of site-specific species within genus Streptomyces can offer insight on the adaptive capacity of strains to inhibit locally coexisting resource competitors within and among these distinct locations.Distribution of bioactive Streptomyces species at different sediment depthsWe further investigated the distribution of antibiotic-producing Streptomyces strains along with the 110-cm sediment depth in different sampling sites. From the sediments that were partitioned according to depth with five categories at 25-cm increments, heterogeneous distributions of bioactive species were observed at deeper sediment with different dominant species in each depth (Fig. 6c). Although, S. parvulus, S. rochei, and S. enissocaesilis strains were ubiquitous in all depths, several species thrive abundantly in specific depths compared to other species. Streptomyces parvulus was the most dominant species in surface sediment. Meanwhile, Streptomyces rochei and S. enissocaesilis strains were more adapted in surface and sub-bottom sediments, respectively. Depth-specific Streptomyces strains were also identified as follows: Streptomyces sp. strain DSD3025 was isolated in subsurface sediments; S. mutabilis strains were abundant in the middle sediment layer; Streptomyces sp. strain DSD1006 and S. pseudogriseolus strain were recovered from sub-bottom sediments; and Streptomyces strain sp. DSD742, S. albus, S. sedi and S. xiamenensis strains were obtained from bottom sediments. High species diversity was positively correlated with increasing sediment depth, where surface sediments are known to be more prone to dispersal and wash-offs by environmental factors such as deep ocean currents42. Furthermore, the depth-specific species identified largely influenced the species richness in varying sediment depth. More

Study regions and farmsTen European regions from boreal to Mediterranean were selected (Supplementary Table 1). They represented major agricultural land uses such as arable crops including horticulture, mixed farming, grassland and perennial crops (vineyards and olives). Within each region, a pool of ~20–40 farms was selected from which 12–20 farms were randomly selected (169 in total) that belonged to the same farm type, produced under homogeneous climatic and environmental circumstances and fulfilled specific criteria regarding their main production branch. In case the selected farms were not willing to participate, we asked other farms from the pool till the sufficient number has been reached. The selected organic farms had all been certified for at least five years. Farmers were asked if they were willing to participate in the study. If they refused, additional random sampling was conducted. In the region NL, 11 organic farms agreed to participate but only three non-organic farms, whereas seven organic farms and 11 non-organic farms were available in the region HU. During the study, one non-organic farmer in the region CH ceased participation.Habitat maps and farm interviewsThe complete area of all selected farms was mapped, using the BioHab method36. Excluded from the farm area were woody and aquatic habitats larger than 800 m2 and summer pastures. Within the farm area, areal and linear habitats were recorded. For an areal habitat, the minimal mapping unit was 400 m2 with a width of at least 5 m. More narrow habitats, between 0.5 and 5 m wide and at least 30 m long, were mapped as linear habitats. Habitats were distinguished in habitat types according to Raunkiær life forms, environmental conditions and management evidence28. Further, a farmland class was assigned to each habitat that described whether the habitat was managed for agricultural production or other objectives such as e.g. nature conservation. In face-to-face interviews following a standardized questionnaire, farmers provided detailed information on field management and yield.Categorization as production fields and semi-natural habitatsBased on the habitat maps and available information about management intensity, we categorized all habitats as either semi-natural habitats or production fields. In agricultural landscapes, these two categories are often not clearly distinguishable. There is a gradient from more intensively managed production fields to less intensively used semi-natural habitats. In addition, a categorization at the local scale can be different from an approach at a European scale (29 and see p. 45 of37). Here, we applied the same criteria for all ten study regions.In all cases, we categorized as production fields: arable crops, intensively managed grasslands (following main plant species observed, management evidence and objectives, with fertilization and/or two or more cuts a year), horticultural crops, and vineyards.We categorized as semi-natural habitats: linear habitats, habitats that were managed for nature conservation objectives, habitats where mainly geophytes, helophytes or hydrophytes were growing, grasslands with woody vegetation (shrubs and/or trees), and extensively managed grasslands (no fertilization, no or one cut a year).Species samplingVascular plant, earthworm, spider and bee species were sampled in all different habitat types of a farm. One plot per habitat type was randomly selected per farm for species sampling. This resulted in 1402 selected habitat plots on 169 farms (Supplementary Table 2). In the selected habitats, species were sampled during one growing season, using standardized protocols19,38. Plant species were identified in squares of 10 × 10 m in areal habitats and in rectangular strips of 1 × 10 m in linear habitats. Earthworms were collected at three random locations of 30 × 30 cm per habitat. First, a solution of allyl isothiocyanate (AITC) was poured out to extract earthworms from the soil. Afterwards, a 20-cm-deep soil core from the same location was hand sorted to find additional specimens. Identification took place in the lab. Spiders were sampled on three dates at five random locations per habitat within a circle of 0.1 m2. Using a modified vacuum shredder, spiders were taken from the soil surface, transferred to a cool box, frozen, or put in ethanol, sorted and identified in the lab. Bees (wild bees and bumble bees) were sampled on three dates, during dry, sunny and warm weather conditions. They were captured with an entomological aerial net along a 100 m long and 2 m wide transect, transferred to a killing jar and identified in the lab.Grouping of species dataSpecies data were pooled per taxa, habitat and region, and three sub-communities were formed with species (1) exclusively found in semi-natural habitats, i.e. unique to semi-natural habitats, (2) exclusively found in production fields, i.e. unique to production fields, and (3) found in both habitat categories i.e. shared by production fields and semi-natural habitats. For calculations of effects over all four taxa, species richness was the sum of the individual taxa species richnesses.Estimating species richnessSpecies richness was estimated using coverage- and sample-size-based rarefaction and extrapolation curves31,39,40. Rarefaction and extrapolation, including confidence intervals (bootstrap method) and sampling coverage, were calculated in R 3.4.041 using package iNEXT42. Detailed information is provided below for each topic.Estimating richness of unique species to compare semi-natural habitats and production fieldsTo legitimately compare the richness of species unique to semi-natural habitats and to production fields, we used the coverage-based method, i.e. we standardized the samples by their completeness30. The point of comparison was determined by the so-called ‘base coverage’ identified by the following procedure31: (1) select the maximum sample coverage at reference sample size (number of sampling units) of the sub-communities under comparison, (2) select the minimum sample coverage at twice the reference sample size of the sub-communities under comparison, (3) identify the maximum of the results from step (1) and step (2) as ‘base coverage’. The species richness estimates were then read off from the species sample-size-based rarefaction and extrapolation curves at the ‘base coverage’ for each sub-community being compared. If zero or exactly one species was unique to a sub-community at the reference sample size, no sample coverage could be calculated. In this case, we set the species richness at 0 or 1, respectively. The species richness estimate of the other sub-community under comparison was then read off at twice the reference sample size on the curve.The ‘base coverage’ was individually defined for each region and each taxonomic group since the mixed effects models used to analyze the data took into account the variation among regions and taxonomic groups.Differences in species richness unique to semi-natural habitats and production fieldsThe difference between the species richness unique to semi-natural habitats and unique to production fields was tested with mixed effects models using package lme4 (Version 1.1-12) in R43. The data were (Sij | β, b, x) ~ Poisson(µij) from i = 1, …, 10 regions. The model is:$${{{rm{ln}}}}left({mu }_{{ij}}right)={beta }_{0}+{beta }_{1}{x}_{1i}+{b}_{1i}$$
(1)
$${b}_{1} sim N(0,sigma 2)$$where ({beta }_{0}) is a fixed intercept, ({beta }_{1}) a fixed effect sub-community ({x}_{1{ij}}) (species unique to semi-natural habitats versus species unique to production fields), b1i are random intercepts for region i. Random effects are normally distributed with mean 0 and variance σ2. The significance of term ({beta }_{1}) was calculated by log-likelihood ratio tests with one degree of freedom. For the models over all four taxa, an additional random intercept was included, i.e. b2j with mean 0 and variance σ2 for j = 1, …, 4 taxa (Fig. 1b).Differences in species richness between organic and non-organic systemsThe comparison between organic and non-organic systems of species unique to semi-natural habitats and to production fields, and of species shared by the two habitat categories, relied on coverage-based extrapolation as described above. Differences between management systems were tested for significance using mixed-effects models with management system ({beta }_{1}) ({x}_{1{ij}}) as fixed effect in (1).Estimating species loss due to conversion of semi-natural habitats to production fieldsTo predict the species loss due to conversion of semi-natural habitats to production fields, we relied on sample-size-based extrapolations31 with species incidence frequencies. We estimated the richness of the species pool for the total number of mapped habitats including the extrapolated species richness unique to semi-natural habitats and unique to production fields, and the observed richness of shared species for each of the four taxa. This species pool provided the basis for the calculation of the species loss or gain (Table 1 and Supplementary Table 7). To model the species richness decrease for any amount of semi-natural habitats converted to production fields, we calculated and drew backward the curve composed of the accumulation curve for species unique to semi-natural habitats, to which the estimated total species richness unique to production fields (constant) and the corresponding gain of species unique to production fields (increases with increasing area of production fields as semi-natural habitats are converted), and the richness of observed shared species (constant) were added. This is the species decrease curve (Supplementary Fig. 2). If started at the observed species richness, this curve corresponds exactly to a species richness curve calculated by a cumulative random removal of semi-natural habitats one by one from the pool of all habitats. The four taxa decrease curves were added for the curve in Fig. 2. Confidence intervals (CI, 95%) shown in Figs. 2 and 3 are calculated by bootstrapping within the calculation of the species accumulation curves (iNEXT42), upper and lower bounds of the 95% CI of the four taxa being added. From the species decrease curve, we read off the predicted species richness for a conversion of 50% and 90% of the semi-natural habitats, and a conversion required to increase production by 10%.As species were sampled in 20% of all mapped habitats on average per region (min. 8%, max. 35%), extrapolated species accumulation curves used to build the species decrease curve were calculated for more than two to three times the reference sample size, which is the suggested range for reliable extrapolation of the species richness estimator31,44. Obviously, the confidence intervals (CI) of the species richness extrapolations here became wide (Supplementary Fig. 4). As we still wanted to show the impact of a conversion of the whole semi-natural area into production fields on the production gain in the ten regions, we used the uncertainty (upper and lower bounds of the 95% CI of the four taxa added) to define two situations in addition to the average case to predict species richness for a 50% and a 90% semi-natural habitat conversion, and a conversion required to increase production by 10%: (1) a worst case situation with the upper bound of the CI of the expected species richness unique to semi-natural habitats, the lower bound of the CI of the expected species richness unique to production fields, and shared species assumed not to be able to survive without semi-natural habitats and considered like species unique to semi-natural habitats (i.e. upper bound); and (2) a best case situation with the lower bound of the CI of the expected species richness unique to semi-natural habitats, the upper bound of the CI of the expected species richness unique to production fields, and the lower bound of the CI of the expected shared species richness.Estimating production gainFarmer interviews delivered an average yield per crop type per farm for the years 2008–2010 (Supplementary Data45 shows details for organic and non-organic systems separately). Farmers indicated yield in kilograms or tons per hectare. This was transformed into energy units, i.e. mega joules per hectare (MJ ha−1) using standard values46. From this, for each region, the average yield (MJ ha−1) was calculated by first multiplying individual crop type yields by the corresponding crop type areas to obtain the production per crop type, then summing up the production of all crop types, and finally dividing this sum by the total area of the crop types. For livestock farms, the fodder production of grasslands was estimated based on the average requirements per livestock unit, accounting for the amount of feed grain, legumes, silage maize and of imported feedstuff. All yields relate to plant biomass production and do not comprise livestock products. The average yield takes into account the relative cover of the different crop types in the regions. Therefore, the conversion of the semi-natural area to production fields was region-specific. The production of certain semi-natural habitats as e.g. olive groves in Spain was not part of the production calculation. The reason is that data on production for semi-natural habitats were mainly not available and/or negligible, e.g. extensively used grassland in CH or in HU, and we decided to apply the same treatment to all the regions. Consequently, in case of olive groves in Spain the effective increase in production is overestimated. To calculate the production gain per region, the production field area added by the conversion of semi-natural habitat area was multiplied by the average yield. In practice, in many regions it may be impossible to convert semi-natural habitat to productive land due to geomorphological constraints and poor soils, and even if land were converted, yields would be much lower than these averages. The results presented here, especially the 90% scenario, are therefore over-optimistic. On the other hand, our calculations are based on the area of semi-natural habitat available for conversion on existing farms, but in some regions other sources of semi-natural land may be available for conversion, e.g. former agricultural land that has been abandoned.Species loss and production gain for three scenariosWe calculated the change of species richness and the production gain under current day production efficiency for two scenarios: (1) a conversion of 90% of the semi-natural area into production fields. The 10% of semi-natural area remaining is considered unsuitable for agricultural use or even impossible to cultivate; (2) a conversion of 50% of the semi-natural area into production fields, and (3) a necessary conversion of the semi-natural area into production fields to achieve a 10% production increase per region.Standardization for organic and non-organic systemsAlthough the overall mapped area, the number of semi-natural habitats, the number of production fields and the average habitat size did not significantly differ between the two management systems (Supplementary Table 5), we standardized the number and size of habitats to the average across both systems per region to compare the species loss and production gain at current day production efficiency in the organic and non-organic systems. The total production in organic and non-organic systems per region was calculated based on the respective yield and the average mapped area of the production fields across both systems as described in section “Estimation of production gain”. The impact on biodiversity was analyzed for the scenario that organic systems should achieve the same level of production as non-organic systems by converting semi-natural habitats to production fields. We calculated the amount of the required area to be converted into production fields and the corresponding species change.Differences between management systems were again tested for significance using mixed-effects models with management system ({{{{rm{beta }}}}}_{1}) ({{{{rm{x}}}}}_{1{{{rm{ij}}}}}) as fixed effect in (1).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Africa’s population is rapidly growing, with its share of the global population projected to increase from 17% in 2020 to 39% by 2100 (ref. 8). The continent is already grappling with low agricultural productivity and food security challenges. Tremendous efforts are needed to increase food production; however, arable land continues to undergo widespread degradation due to issues such as nutrient mining, erosion, overgrazing and pollution. Climate change and more frequent weather extremes, such as floods and droughts, further degrade land and reduce agricultural productivity.Some efforts to counteract low productivity, however, can increase greenhouse gas emissions and derail efforts to meet global climate targets. Poor water management, fertilizer application and residue burning in rice production are, for example, major sources of potent greenhouse gases such as methane and nitrous oxide9,10. To ensure that the United Nations sustainable development goals and the African Union’s Agenda 2063 for food and water security are realized at minimal environmental cost, science-based land management practices are needed to decouple agricultural productivity from greenhouse gas emissions.
Credit: majimazuri21/PixabayThe Agriculture, Forestry and Other Land Uses (AFOLU) sector contributes the largest share of greenhouse gas emissions in Africa11. Thus, developing large-scale agronomic, livestock and forest management practices that increase productivity and reduce emissions is key to achieving enhanced production and environmental sustainability. However, it is impossible to effectively manage greenhouse gas emissions if there is limited capacity to quantify them in Africa.Improved data infrastructure and research are needed to quantify emissions associated with specific land management practices under different land uses. Similarly, land use mitigation strategies should be informed by existing and potential future land use changes and their impact on greenhouse gas emissions under different climate scenarios. However, past studies that examined land use changes at various temporal scales mainly used coarse resolution satellite imagery and suffered from limited availability or poor-quality of data, partly due to cost. Such challenges have resulted in limited knowledge of land management practices that reduce greenhouse gas emissions while increasing agricultural productivity.Improved greenhouse gas observation networks and in situ measurements12 will enable the development of country-specific emission factors (IPCC tier 2/3)13 and quantification and management of land use specific greenhouse emissions. It will reduce uncertainties in emissions inventory data on Agriculture, Forestry and Other Land Uses14, which are currently estimated using emission factors extracted from default value databases (tier 1 methodologies).Free earth observation data, such as those from the European Space Agency and United States Geological Surveys, are becoming increasingly available. Together with improvements in cloud-based computing infrastructure, this presents an opportunity to advance research into current and future land use and vegetation dynamics. Coupled with accurately quantified greenhouse gas emissions, this can support current and future land management practices that contribute to mitigation and adaptation objectives of countries. More

Sample collection of Liolophura japonica and the genetic diversity of COI barcoding regionTo examine genetic lineage divergence within L. japonica on the northwestern Pacific coast, we newly collected a total of 342 L. japonica samples from 12 sampling localities in the intertidal coasts of the Korean Peninsula and Japanese Archipelago (Fig. 1; Table S1). From the collected L. japonica samples, we amplified the COI barcoding region using PCR, and then sequenced the 635-bp PCR products. As a result, a total of 75 COI haplotypes based on COI sequences obtained from 342 individuals of L. japonica were detected via the present study (Table S2). In addition to this, we extracted 31 COI haplotypes based on COI sequences from 127 individuals of L. japonica (also known as Acanthopleura japonica) previously reported in the NCBI GenBank database, consisting of two Japanese and 29 Chinese COI haplotypes. Finally, we gathered 106 COI haplotypes from 469 L. japonica individuals collected in 15 localities of South Korea, Japan, and southern China (Tables S1, S2). The average haplotype (h) and nucleotide diversities (π) were 0.808 and 0.04936, respectively; the highest haplotype diversity was observed in Tsushima (TS; h = 0.963), and the highest nucleotide diversity was found in Wando (WD; π = 0.04581), located in the South Sea of the Korean Peninsula. As shown in Table S3, the population distribution pattern of COI haplotypes revealed that all collection sites had site-specific haplotype(s) except for Busan (BS), Wando (WD), Sinan (SA), and Jeju Island (JJ). The most abundant haplotype was A1, which was found in 170 (39.4%) out of the COI sequences obtained from 469 L. japonica individuals.Figure 1A map showing sampling localities and photos of a habitat landscape and wild samples of Liolophura japonica inhabiting coastal areas of the Korean Peninsula (N = 249), the Japanese Archipelago (N = 57), and southern China (N = 125) in the northwestern Pacific Ocean. (a) A map showing twelve direct sampling localities for L. japonica in coastal areas of the northwestern Pacific Ocean. The sampling localities of one southern Chinese (ZJ) and two Japanese (EH and MY) previously catalogued haplotype sequencing studies retrieved from NCBI are also depicted. Table S1 and S2 contain more accurate information on the populations and individuals. The basic map is from a free map providing site (https://d-maps.com), which is modified with Adobe Illustrator v.25.2. (https://www.adobe.com). (b) Photos of a habitat landscape and wild samples of L. japonica, taken from Seogwipo-si, Jeju Island, South Korea, photographed by Mi Young Yeo, Bia Park, and Cho Rong Shin. The photos were edited using Adobe Photoshop v.22.2 (https://www.adobe.com).Full size imagePhylogenetic and population genetic analyses based on COI
We constructed a nucleotide sequence alignment set with 106 COI haplotypes of L. japonica (Data S1), and identified 95 polymorphic sites (15.0%, Table S4) and 68 parsimoniously informative sites (10.7%). To elucidate phylogenetic relationships among the populations of L. japonica, we performed molecular phylogenetic analyses, including maximum likelihood (ML), Bayesian inference (BI), and neighbor-joining (NJ) analyses, based on these 106 COI haplotypes with the outgroup Acanthopleura spinosa (Fig. 2a, Figs. S1, S2). The resultant phylogenetic trees clearly revealed the existence of three distinct genetic lineages within the monophyletic group of L. japonica (100 BP in ML, 1.00 BPP in BI, and 100 BP in NJ): Lineage N (91 BP, 1.00 BPP, and 100 BP), Lineage S1 (79 BP, 0.82 BPP, and 98 BP), and Lineage S2 (98 BP, 1.00 BPP, and 100 BP). Among these three genetic lineages, Lineages S1 and S2 were grouped with high node confidence values (94 BP, 1.00 BPP, and 95 BP). We additionally conducted a phylogenetic network analysis using a neighbor net algorithm without an outgroup (Fig. 2b), which confirmed that these sequences were distinctly divided into three genetic lineages, in agreement with the topology of the rooted phylogenetic trees (Fig. 2a, Fig. S1).Figure 2Phylogenetic, TCS network, and PCoA analyses based on 106 COI haplotypes from 469 individuals of Liolophura japonica inhabiting coastal areas of the northwestern Pacific Ocean, suggesting the existence of the three different genetic lineages: Lineage N, Lineage S1, and Lineage S2. (a) Maximum likelihood tree showing the three different genetic lineages for L. japonica: Lineage N members are most likely from the populations inhabiting a wide range of South Korea and Japan, Lineage S1 members from the populations inhabiting southern coastal areas of South Korea and Japan only, and Lineage S2 members from the southern Chinese population. As shown in Fig. S1, Acanthopleura spinosa was used as an outgroup. Numbers on branches indicate node confidence values: BP in ML, BPP in BI, and BP in NJ in order. (b) A phylogenetic network reconstructed using the neighbor net algorithm without an outgroup, showing three different genetic lineages for L. japonica inhabiting the northwestern Pacific coast: Lineages N, S1, and S2. The COI sequence alignment set used is shown in Data S1. Detailed information of the 106 COI haplotypes used in this phylogenetic analysis is summarized in Table S1 and S2. (c) An unrooted TCS network showing three distinct genetic clusters, corresponding to Lineages N, S1, and S2. Three different genetic groups correspond to the three genetic lineages shown in the phylogenetic tree (a), respectively. The haplotype frequency is displayed by the circle size. (d) A two-dimensional PCoA plot showing the three distinct genetic groups corresponding to Lineages N, S1, and S2 shown in the phylogenetic tree (a). The score on the first two axes (Axis 1 = 79.05% and Axis 2 = 15.32%) from the matrix of genetic distances estimated with the 106 COI haplotypes are indicated.Full size imageConsistently, the TCS network analysis (Fig. 2c) and principal coordinate analysis (PCoA) (Fig. 2d) showed the existence of three distinguished genetic groups among L. japonica, in accordance with the phylogenetic analyses (Fig. 2a − b). The TCS network (Fig. 2c) revealed that Lineages S1 and S2 were separated by 18 mutation steps, which is far shorter than the distance between Lineages N and S1 (37 mutation steps) or between Lineages N and S2 (60 mutation steps), indicating that Lineages S1 and S2 have a close affinity and had only recently diverged from each other. The overwhelming dominancy of the A1 haplotype implies a recent and rapid population expansion of Lineage N. In addition to A1, it was found that haplotypes B2 for Lineage S1 and C21 for Lineage S2 were dominant. In the PCoA plot (Fig. 2d), the three genetic groups of L. japonica were also observed, as in the phylogenetic (Fig. 2a,b) and TCS network (Fig. 2c) analyses. Lineage N was distantly located from Lineages S1 and S2, while Lineages S1 and S2 were spatially much closer.Sample collection of L. japonica and the genetic diversity of 16S rRNA
The 342 individuals of L. japonica from 12 localities in the intertidal coasts on the Korean Peninsula and Japanese Archipelago (Fig. 1) were subjected to PCR amplification of a partial region of 16S rRNA (506 bp) (Tables S5, S6). Of these, only 299 samples were successfully amplified and sequenced. Based on 299 individual 16S rRNA sequences, a total of 23 16S rRNA haplotypes of L. japonica were detected (Tables S5, S6). Combined with 11 haplotypes extracted from 16S rRNA sequences of 125 L. japonica individuals known previously in southern China, we totaled 34 16S rRNA haplotypes from 425 L. japonica individuals in 13 collection localities. The average haplotype (h) and nucleotide (π) diversities were 0.702 and 0.02093, respectively; the highest haplotype diversity was found in Geojedo (GJ; h = 0.833), and the highest nucleotide diversity in Wando (WD; π = 0.02244), located in the South Sea of the Korean Peninsula. Overall, the average haplotype and nucleotide diversities of 16S rRNA were lower than those of COI (Table S1). As shown in Table S7, the population distribution pattern of 16S rRNA haplotypes revealed that most of the collection sites had site-specific haplotype(s), except for BS, GJ, WD, SA, JJ, and TT. The most abundant haplotype was RA1, which was found in 186 (48.1%) out of the 16S rRNA sequences obtained from 425 L. japonica individuals.Phylogenetic and population genetic analyses based on 16S rRNA
We constructed a nucleotide sequence alignment set with 34 16S rRNA haplotypes of L. japonica (Data S2), and identified 35 polymorphic sites (6.9%; Table S8) and 24 parsimoniously informative sites (4.7%). Phylogenetic analyses, including ML, BI, and NJ analyses, were conducted with the outgroup Acanthopleura echinata (Table S6). The resultant phylogenetic trees (Fig. S2) and unrooted phylogenetic network (Fig. 3a) consistently supported the three distinct genetic lineages of L. japonica, with the phylogenetic relationship between Lineages S1 and S2 being much closer than those inferred from the results of COI (Fig. 2a,b; Fig. S1). The TCS network (Fig. 3b) revealed that Lineages S1 and S2 were closely connected with only with 4–5 mutation steps between them, while Lineages N and S1 or Lineages N and S2 were remotely distanced by 18 mutation steps. Also, the overwhelming dominance of the RA1 haplotype implied a recent and rapid population expansion of Lineage N. In addition to RA1, haplotypes RB1 for Lineage S1 and RC1 and RC2 for Lineage S2 were dominant (Fig. 3b; Table S7). Consistent with this, in the PCoA plot (Fig. 3c), the three genetic groups of L. japonica were spatially separated. Lineage N was distantly located apart from Lineages S1 and S2, while Lineages S1 and S2 were spatially much closer.Figure 3The results of phylogenetic and population genetic analyses based on 34 16S rRNA haplotypes from 425 individuals of Liolophura japonica inhabiting coastal areas of the northwestern Pacific Ocean. (a) Phylogenetic network reconstructed using the neighbor net algorithm, showing three different genetic lineages for L. japonica: Lineage N, Lineage S1, and Lineage S2. The 16S rRNA sequence alignment set used is shown in Data S2. Detailed information of 34 16S rRNA haplotypes used in these analyses is summarized in Table S5 and S6. (b) An unrooted TCS network. There are distinctly observed three different genetic groups, corresponding to the three genetic lineages shown in the phylogenetic network (a). The haplotype frequency is displayed by the circle size. (c) A two-dimensional PCoA plot showing the three distinct genetic groups, corresponding to Lineage N, Lineage S1, and Lineage S2. The score on the first two axes (Axis 1 = 87.77% and Axis 2 = 4.4%) from the matrix of genetic distances estimated with the 34 16S rRNA haplotypes are indicated.Full size imageExamination of species discrimination of L. japonica based on COI and 16S rRNA
Using the Automatic Barcode Gap Discovery (ABGD), we performed distribution of pairwise genetic divergences, ranked pairwise difference, and automatic partition analyses based on COI and 16S rRNA of L. japonica, respectively (Fig. 4a–c), which confirmed that there were distinct barcoding gaps between intraspecific and interspecific variations, strongly supporting the possibility of species discrimination of L. japonica. the COI-based analysis yielded two different barcoding gaps, while the 16S rRNA-based analysis revealed only a single barcoding gap (Fig. 4a–c). The results of automatic partition at each value of the prior intraspecific divergence (P) divided L. japonica into three groups by COI and two groups by 16S rRNA, respectively (Fig. 4a–c). We also implemented two DNA taxonomy approaches to evaluate the possibility of species discrimination based on COI: the general mixed Yule coalescent (GMYC) approach (Fig. S3) and a Bayesian implementation of a Poisson Tree Processes model (bPTP) (Fig. S4). The results consistently and robustly supported the possibility that L. japonica can be divided into three different species, as shown in the results of ABDG (Fig. 4a–c).Figure 4Distribution of pairwise genetic divergences, ranked pairwise difference, and automatic partition based on COI and 16S rRNA haplotypes of Liolophura japonica and a COI-based NJ tree showing the phylogenetic relationship with a congeneric species L. tenuispinosa. (a) Distribution patterns of pairwise genetic divergences observed in COI and 16S rRNA for L. japonica. The horizontal axis represents intervals of pairwise Kimura-2-parameter (K2P) genetic distance in percentage, and the vertical axis represents the number of individuals associated with each distance interval. (b) The results of ranked pairwise differences based on COI and 16S rRNA, ranked by ordered value, which is similar to the distribution of pairwise genetic divergence in (a). The horizontal axis indicates a ranked ordered value based on K2P genetic distance, and the vertical axis represents the K2P genetic distance in percentage. (c) The results of automatic partition analyses based on COI and 16S rRNA. The horizontal axis represents the prior maximum intraspecific divergence (P), and the vertical axis represents the number of groups inside the partitions (primary and recursive). (d) A COI-based NJ tree with L. tenuispinosa. Refer to Fig. S3 and Data S3.Full size imageThe molecular variance analyses using analysis of molecular variance (AMOVA), based on COI and 16S rRNA, were conducted to evaluate the degree of genetic differentiation among Lineages N, S1, and S2 (Tables S9, S10). According to the results, supposing that there are three genetic lineages (N, S1, and S2) or two genetic lineages (N and S1/S2), almost all variation in both cases is attributed to variation among groups (= among lineages), whereas variations within populations (within lineages) exhibit negative values in common. We confirmed that there was a high degree of genetic differentiation among Lineages N, S1, and S2, which supports the results of the COI barcoding gap analysis shown in Fig. 4a–c, although this was not statistically significant (P  > 0.05; Table S9). When we assumed only two genetic groups, Lineages N and S1/S2, the genetic differentiation between the two groups was statistically significant (P  0.05) (Tables S9 and S10). The discrepancy between the number of barcoding gaps inferred from COI and 16S rRNA may have been affected by different gene evolutionary rates of the molecular markers11; nucleotide substitution rate of 16S rRNA is known to be generally slower than that of COI (which is especially fast in the third codon positions: 105 out of 127 polymorphic sites). When an ML tree was constructed based on 22 polymorphic sites, which are found only in the first and second codon positions of COI that are much more conserved than the third codon position, the three genetic lineages were retained in the resultant tree (Fig. S5), but Lineage S2 was nested within Lineage S1, as in the trees inferred from 16S rRNA (Fig. S2). Reflecting the powerful resolution of the COI barcoding marker well known from animals12 and the high degree of variation among the three genetic lineages (Fig. 4, Figs. S3, S4), we suggested that L. japonica could be categorized into three different species: L. koreana, Yeo and Hwang, sp. nov. for Lineage N, L. japonica for Lineage S1, and L. sinensis Choi, Park, and Hwang, sp. nov. for Lineage S2. To examine whether it is reasonable to give these a species-level taxonomic status, as shown in Fig. 4d, we reconstructed a COI-based NJ tree with one congeneric species L. tenuispinosa13, which was originally described as a subspecies-level taxon of L. japonica14,15 and was then revised as an independent species closely related to L. japonica by Saito & Yoshioka16 in 1993. The resultant tree (Fig. S6 and Data S3) showed that L. tenuispinosa forms a sister group with L. japonica (Lineage S1). This likely indicates that L. koreana and L. sinensis have taxonomic status as independent species.Morphological comparison and geographical distribution of the three Liolophura speciesWe compared morphological characteristics among Liolophura koreana, sp. nov. (Lineage N), L. japonica (Lineage S1), and L. sinensis, sp. nov. (Lineage S2). Their morphological appearances are shown in Fig. 5a–c, which indicated that black spots on the tegmentum (Fig. 5d–e) and shapes of spicules on the perinotum (Figs. 5f–k, 6e–f) represent key morphological characteristics to distinguish them from each other. Although black dots in pleural areas, which are between the middle and lateral areas of the tegmentum on valves II–VII (or VIII), are commonly shared in all three lineages (Fig. 5a–c), other black spots on the valves exhibit a high degree of variation in morphology (Fig. 5a–c, Fig. S7). Herein, we described a new species of genus Liolophura, that is, L. koreana Yeo and Hwang from South Korea and Japan, with detailed descriptions of morphological characteristics observed by light microscopy (M205, Leica Camera AG, Germany) and FE-SEM (SU8220, Hitachi, Japan). In addition, we suggested the divergence of a new species, L. sinensis Choi, Park, and Hwang from southern China, with simple remarks based on distinct genetic difference (mainly COI barcoding gaps), with possible unique morphological characteristics as follows.Figure 5Morphological comparison of Liolophura koreana, sp. nov., L. japonica, and L. sinensis, sp. nov. (a–c) Photos of dorsal views of the individuals belonging to L. koreana (Lineage N), L. japonica (Lineage S1), and L. sinensis (Lineage S2) in order. (d,e) Morphological comparison of pleural and lateral black spots on valves III and IV of the tegmentum of L. koreana (d; holotype) and L. japonica (e). (f,g) Morphological comparison of spicules on the perinotum of L. koreana (f; holotype) and L. japonica (g). (h–k) Morphological comparison of the spicule of L. koreana (h,i; paratype) and L. japonica (j,k) in lateral and dorsal views. The scale bar marks 2.0 mm (d,e), 1.0 mm (f,g), and 0.5 mm (h–k). The photos were edited using Adobe Photoshop v.22.2 (https://www.adobe.com).Full size imageFigure 6Microstructural comparison of Liolophura koreana, sp. nov. and L. japonica using field emission scanning electron microscopy (FE-SEM). (a,b) Middle and lateral areas on the tegmentum of the holotype of L. koreana. (c,d) Middle and lateral areas on the tegmentum of L. japonica. Arrows indicate that morphological difference of the posterior valve margin of the valve II between two species. The scale bar marks 1.0 mm. (e,f) The occurrence frequency, and shape and structure differences of the spicules on the perinotum between the holotype of L. koreana (e) and L. japonica (f). The scale bar marks 1.0 mm and 0.2 mm, respectively. The photos were edited using Adobe Photoshop v.22.2 (https://www.adobe.com).Full size image
Liolophura koreana Yeo and Hwang, sp. nov. (Figs. 5, 6; Figs. S7, S8)(urn:lsid:zoobank.org:act:4418355E-F55C-44FA-B4CE-585589FDCD23).Type specimens examined[Holotype] SOUTH KOREA: 1 specimen, Jeju-do, Seogwipo-si, Seongsan-eup, Seopjikoji, 3.XI.2020, UW Hwang, B Park & CR Shin (LEGOM040501); [Paratypes] SOUTH KOREA: 1 specimen, Gyeongsangbuk-do, Pohang-si, Guryongpo-eup, Janggil-ri, 27.VII.2008, UW Hwang (LEGOM040502); 3 specimens, Gyeongsangbuk-do, Ulleung-gun, Seo-myeon, Namyang-ri, Ulleungdo Island, Namtong tunnel, 12.VI.2007, UW Hwang (LEGOM040503–0505); 1 specimen, Gyeongsangbuk-do, Ulleung-gun, Namyang-ri, Ulleungdo Island, Namyang tunnel, 5.X.2007, UW Hwang (LEGOM040506); 2 specimens, Gyeongsangnam-do, Geoje-si, Nambu-myeon, Dapo-ri, 28.IV.2009, MY Yeo (LEGOM040507,0508); 4 specimens, same data as the holotype (LEGOM040509–0512); 2 specimens, Jeollanam-do, Yeosu-si, Hwajeong-myeon, Sado-ri, Sado Island, 8.IV.2008, MY Yeo (LEGOM040513,0514); 4 specimens, same data as the holotype (LEGOM040515–0518); JAPAN: 6 specimens, Tottori Prefecture, Hakuto, 24.V.2009, UW Hwang (LEGOM040519–0524); 1 specimen, Tottori Prefecture, Iwato, 25.V.2009, UW Hwang (LEGOM040525).DescriptionBody small-sized and broad oval- to oval-shaped (Fig. 5a; Fig. S7); length 3.9 (1.9–12.3) mm and width 2.4 (1.2–7.1) mm. Tegmentum entirely brown (dark brown or black entirely, or each valve with black line anteriorly or white line laterally), with black dots on the pleural areas of valves II–VII (or VIII) (Fig. 5a,d, Fig. S7); articulamentum entirely black (dark brown); whitish and blackish spicules on the perinotum scattered irregularly, sometimes forming a band besides each valve (Fig. 5a, Fig. S7). Surface of the tegmentum in middle and lateral areas as in Fig. 6a,b and Fig. S8; posterior margin of the head valve nearly straight; dorsal shape of intermediate valves round-backed and side slopes slightly convex; the posterior valve margin with a distinct central apex, its shape subtriangular to triangular (rounded or linear), particularly valve II, mainly with a strong projection (Fig. 6a). Perinotum covered with large, solid, slightly curved, and obtusely pointed spicules (rarely with smooth and radial ribbed spicules apically), its density relatively lower than that of L. japonica (Figs. 5f,h,i, 6e).DistributionSouth Korea, Japan; avobe 33°24′ N (Seogwipo, JJ) in South Korea and TT and MY in Japan (Fig. 7).Figure 7Geographical distribution of Liolophura koreana, sp. nov., L. japonica, and L. sinensis, sp. nov. inhabiting coastal areas of the northwestern Pacific Ocean. A COI-based map showing geographical distribution of L. koreana, L. japonica, and L. sinensis on the northwestern Pacific coast. L. koreana are found in a wide range of South Korea and Japan above ca. 33°24′ N (JJ), L. japonica in mainly southern coastal areas of South Korea and Japan below ca. 35°53′ N (TT), and L. sinensis in ZJ of southern China around ca. 27°02′ N–28°00′ N. The sympatric distribution of L. koreana and L. japonica is shown between 33°24′ and 35°53′ N. Table S1–S3 contain the full names of localities and detailed haplotype information. The question mark indicates that collection of Liolophura samples from such coastal areas in Japan is required to clarify distribution patterns of L. koreana and L. japonica in the East Sea (= Sea of Japan). The basic map was obtained from a free map-providing site (https://d-maps.com), which was modified using Adobe Illustrator v.25.2. (https://www.adobe.com).Full size imageHabitatThis new species appears to be attached to rocks in coastal areas with strong waves, or a calm inner shore in the northwestern Pacific Ocean (Fig. S9).EtymologyThe species is named per the locality of the new species.RemarksWe found that Liolophura koreana, sp. nov. has no black spots on lateral areas of the tegmentum (Fig. 5d), and large, slightly curved, and obtusely pointed spicules on the perinotum compared to those of L. japonica (Figs. 5f,h,i, 6e). On the other hand, L. japonica from southern South Korea and southern Japan has two black spots on the lateral areas of valves II–VII (or VIII) (Fig. 5e), and small, almost straight, and cylindrical spicules compared to those of L. koreana (Figs. 5g,j,k, 6f).As shown in Fig. 7, L. koreana (Lineage N) was observed in all the South Korean and Japanese populations examined here, except for the EH population in Japan (refer to Tables S1, S2), which were found from JJ at 33°24′ N to MY at 38°32′ N. On the other hand, L. japonica (Lineage S1) was found from the southern coastal areas of South Korea and Japan, which were found only between JJ at 33°24′ N and TT at 35°53′ N. Interestingly, we found only L. koreana north from the latitude of 35°10′ N (BS) such as UL/DD (37°24′ N) and PH (36°02′ N) in South Korea. In Japan, there was found only L. koreana at MY (38°32′ N) too, but it remains to be explored to clarify its distribution range in Japan with much more sample collections covering northern Japanese coastal areas through further study. It was also confirmed that L. koreana and L. japonica show a sympatric distribution pattern between JJ at 33°24′ N and TT at 35°53′ N on the southern coastal area of South Korea and Japan.
Liolophura sinensis Choi, Park, and Hwang, sp. nov. (Fig. 5c).(urn:lsid:zoobank.org:act:72DF7E75-1853-4F23-AC12-3AB8CD054187).Type specimens examined[Holotype] CHINA: 1 specimen, Zhejiang Province, Dongtou Island, 27°49′57.44″ N, 121°10′19.13″ E, 2017; [Paratypes] CHINA: 1 specimen, Beiji Island, 27°37′08.82″ N, 121°11′47.82″ E, 2017; 1 specimen, Beilongshan Island, 27°40′08.56″ N, 121°58′51.56″ E, 2017; 1 specimen, Chaishi Island, 27°25′40.36″ N, 121°04′54.45″ E, 2017; 1 specimen, Daleishan Island, 27°29′39.48″ N, 121°05′24.50″ E, 2017; 1 specimen, Daqu Island, 27°47′29.92″ N, 121°05′23.97″ E, 2017; 1 specimen, Dazhushi Island, 27°49′12.87″ N, 121°12′48.74″ E, 2017; 1 specimen, Dongce Island, 27°45′32.04″ N, 121°09′01.38″ E, 2017; 1 specimen, Dongxingzai Island, 27°02′40.36″ N, 121°02′47.98″ E, 2017; 1 specimen, Houjishan Island, 27°28′26.96″ N, 121°07′40.71″ E, 2017; 1 specimen, Luxi Island, 27°59′33.43″ N, 121°12′50.70″ E, 2017; 1 specimen, Nanji Island, 27°27′30.57″ N, 121°03′06.28″ E, 2017; 1 specimen, Nanpanshan Island, 28°00′15.29″ N, 121°15′33.62″ E, 2017.DistributionSouthern China; Zhejiang Province around ca. 27°02′–28°00′ N (Fig. 7).HabitatThis new species appears to be attached to rocks in coastal areas of the northwestern Pacific Ocean.EtymologyThe species is named after its locality.RemarksL. sinensis, sp. nov. was examined and established mainly by molecular data such as the COI barcoding gap (Fig. 4a–c) presented in this study and photos provided from Prof. Yong-Pu Zhang (Wenzhou University, Zhejiang Province, China) without any direct real sample observation. The morphology of this new species is very similar to that of a previously known species, L. japonica. L. sinensis from southern China has black spots on the lateral areas of valves II–VII similar to L. japonica, but black bowtie-shaped spots anterior to valves II–VII (Fig. 5c) are a unique characteristic for L. sinensis. L. sinensis (Lineage S2) was found in Zhejiang Province (ZJ) in southern China, around ca. 27°02′–28°00′ N (Fig. 7; Tables S1, S5).Demographic history and divergence time estimation analysesMismatch distribution analyses (MDA) based on COI were performed for L. koreana, L. japonica, and L. sinensis, respectively. The MDA results (Fig. 8a) showed a unimodal curve for each of the three lineages. In addition, when neutrality tests were performed with COI and 16S rRNA (Table S11), all three showed statistically significant negative values in both Tajima’s D and Fu’s Fs, except for the Tajima’s D values in COI (L. sinensis) and 16S rRNA (L. japonensis and L. sinensis), implying that these had experienced population expansions. Bayesian skyline plot (BSP) analyses with COI (Fig. 8b) were performed to examine the fluctuation patterns in effective population sizes for L. koreana, L. japonica, and L. sinensis, respectively. The effective population sizes of L. japonica and L. sinensis had gradually grown between ca. 100 Ka and ca. 50 Ka, while those of L. japonica had grown between ca. 80–50 Ka, L. sinensis had grown between 100–60 Ka, and that of L. koreana had begun to rapidly expand ca. 85 Ka and ceased ca. 75 Ka. This indicated that population expansion had occurred more dramatic in L. koreana than in L. japonica and L. sinensis, following the last interglacial age, called the Eemian (129–116 Ka). As shown in Fig. 8c and Fig. S10, according to the molecular clock analysis by the BEAST program, it was estimated that L. japonica and L. koreana shared their most recent common ancestor about 3.37 Ma, around the mid-Pliocene warm period (3.30–3.00 Ma), before the extensive glaciation in the late Pliocene (ca. 3.00 mya). L. japonica and L. sinensis likely diverged around 1.84 Ma, around the beginning stage of the Early Pleistocene Transition (EPT; 1.85–1.66 mya). The augmentation of haplotype diversity in L. japonica, L. sinensis, and L. koreana might have intensified in the interglacial stages during the late-middle (0.35–0.126 Ma) and late Pleistocene (0.126–0.012 Ma), before the last glacial maximum (LGM: 0.026–0.019 Ma).Figure 8The results of mismatch distribution analyses (MDA), Bayesian skyline plots (BSPs), and molecular clock analysis performed with COI haplotypes for Liolophura koreana, sp. nov., L. japonica, and L. sinensis, sp. nov. (a) MDA plots resulting in a unimodal curve for L. koreana, L. japonica, and L. sinensis. Dotted lines indicate the observed distribution of mismatches, and solid lines represent the expected distribution under a demographic expansion model. (b) BSP results showing the demographic history of population expansions of L. koreana, L. japonica, and L. sinensis. The graph in gray depicts sea level changes during the last 330 Ka. (c) Time-calibrated Bayesian tree reconstructed using BEAST with the inference of ancestral areas under the Bayesian binary MCMC (BBM) model implemented in RASP ver 3.2. Ancestral areas were hypothesized based on the distribution range of the fossil records of Mopalia and the contemporary distribution of L. koreana, L. japonica, and L. sinensis. LGM indicates the last glacial maximum (0.026–0.019 Ma; blue vertical bar) and three interglacial periods are indicated by light green boxes during the late-middle and late Pleistocene. The pictures were edited using Adobe Illustrator v.25.2. (https://www.adobe.com).Full size image More