Philippa Kaur

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

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

We consider public goods games played iteratively over a fixed connected network. The vertices of the network represent the players and the edges represent neighboring connections5,10,11,12. The dynamics evolve in discrete rounds. In each round, each player chooses a behavior that minimizes its cost, where the player’s cost is affected by its own behavior and the behaviors of its neighbors.Our main model includes three behavior types, namely, defection, hypocrisy, and cooperation, in which those who hardly contribute to the social welfare, i.e., defector and hypocritical players, face the risk of being caught and punished by their neighbors who are non-defectors. The level of risk together with the extent of punishment is captured by a notion that we call “social-pressure”. The main result is that adjusting the level of social-pressure employed against hypocritical players compared to the one employed against defectors can have a dramatic impact on the dynamics of the system. Specifically, letting the former level of social-pressure be within a certain range below the latter level, allows the system to quickly transform from being composed almost exclusively of defectors to being fully cooperative. Conversely, setting the level to be either too low or too high locks the system in a degenerate configuration.As mentioned, our main model assumes that non-defectors induce mild social-pressure on the defectors among their neighbors. This implicitly assumes that inducing the corresponding social-pressure is beneficial (e.g., allows for a social-upgrade), although other explanations have also been proposed21. To remove this implicit assumption we also consider a generalized model, called the two-order model, which includes costly punishments. Consistent with previous work on the second-order problem, e.g.,23,25,26,27,36,40, this model distinguishes between first-order cooperation, that corresponds to actions that directly contribute to the social welfare, and second-order cooperation, that corresponds to applying (costly) social-pressure, or punishments, on others. As in the main model, the level of punishment employed against first-order defectors may differ from that employed against second-order defectors. We identify a simple criteria for the emergence of cooperation: For networks with minimal degree (Delta), cooperation emerges when two conditions hold. The first condition states that the cost (alpha _2) of employing punishments against second-order defectors should be smaller than the corresponding punishment (beta _2) itself, i.e., (alpha _2 More

Technical criterion for ecosystem status evaluationOn March 3, 2015, the Ministry of Ecology and Environment of the People’s Republic of China approved the “Technical Criterion for Ecosystem Status Evaluation” as the national environmental protection standard. This standard is based on the former standards released in 2006, and 48 relevant documents from 2006 to 2012 were searched to propose new standards and factor weights based on actual utilization effects and expert guidance. The eco-environment assessment uses a comprehensive index (eco-environmental status index, EI) to reflect the overall state of the regional eco-environment. The indicator system includes the biological abundance index, vegetation coverage index, river density index, land stress index, and pollution loading index. These indexes reflect the abundance of organisms in the evaluated area, the level of vegetation coverage, the abundance of water, the intensity of land stress, and the extent of the pollution load. Each indicator was calculated according to its weight to obtain an eco-environment assessment map (Table 1). All parameters involved in the calculation are derived from this standard.Table 1 Weights of the evaluation indicators.Full size tableThe calculation of the eco-environment status is as follows:$$begin{aligned} EI & = 0.35*biological ; abundance ; index + 0.25*vegetation ; coverage ; index hfill \ & quad + 0.15*river ; density ; index + 0.15*(1 – land ; stress ; index) hfill \ & quad + 0.1*(1 – pollution ; loading ; index) hfill \ end{aligned}$$
(9)
Biological abundance indexThe biological abundance index refers to the number of certain organisms in this area. The calculation method is as follows:$$Biological , abundance , index , = , left( {BI , + , HQ} right)/2$$
(10)
In this formula, BI is the biodiversity index and HQ is the habitat quality index. When the biodiversity index does not have dynamic data updates, the change in the biological abundance index is equal to the change in the HQ.Biodiversity is a general term for the complexity of species and their genetic variation and ecosystems in space over time. Biodiversity plays an important role in maintaining soil fertility, ensuring water quality, regulating the climate, stabilizing the environment, and maintaining ecological balance.The BI method is as follows:$$BI = NPP_{mean} *F_{pre} *F_{tem} *(1 – F_{alt} )$$
(11)
NPPmean is the net primary productivity. Fpre is the annual average precipitation. Ftem is the temperature parameter. Falt is the altitude parameter.NPP refers to the amount of organic matter accumulated per unit area and unit time of green plants. NPP is the remainder of the total amount of organic matter produced by photosynthesis after deducting autotrophic respiration and is usually expressed as dry weight. In this study, the estimation of NPP was based on the absorbed photosynthetically active radiation (APAR) and actual light-use efficiency (LUE) (ε) of the CASA ecosystem model40. The CASA model is a process-based remote sensing model that couples ecosystem productivity and soil carbon and nitrogen fluxes, driven by gridded global climate, radiation, soil, and remote sensing vegetation index datasets41. The model can be expressed generally as follows:$$NPP(x,t) = APAR(x,t)*varepsilon (x,t)$$
(12)
The entire study area is divided into 11,303 pixels on a 30 * 30 m grid. x indicates the location of each pixel, and t indicates time; the data were collected once a month. APAR(x,t) represents the photosynthetically active radiation absorbed by pixel x in that month (gC * m−2* month−1). Ɛ(x, t) is LUE (gC * MJ−1) of the vegetation42.Estimation of the fraction of APAR using RS data is based on the reflection characteristics of the vegetation on the infrared and near-infrared bands. The value of APAR is determined by the effective radiation of the sun and the absorption ratio of the vegetation to the effective photosynthetic radiation. The formula is as follows:$$APAR(x,t) = SOL(x,t)*FPAR(x,t)*0.5$$
(13)

where SOL(x,t) represents the total amount of solar radiation at pixel x in month t, FPAR(x,t) represents the absorption ratio of the vegetation layer to the incident photosynthetically active radiation, and a constant of 0.5 indicates the ratio of the effective solar radiation that can be utilized by the vegetation to the total solar radiation.Since there is a linear relationship between FPAR and NDVI within a certain range, this relationship can be determined according to the maximum and minimum values of a certain vegetation type NDVI and the corresponding FPAR maximum and minimum values.$$FPAR(x,t) = frac{{(NDVI(x,t) – NDVI_{i,min } )}}{{(NDVI_{i,max } – NDVI_{i,min } )}}*(FPAR_{max } – FPAR_{min } ) + FPAR_{min }$$
(14)

where NDVImax and NDVImin correspond to the NDVI maximum and minimum values of the ith planting type, respectively. There is also a good linear relationship between FPAR and the simple ratio index (SR) of vegetation, which is represented by the following formula:$$FPAR(x,t) = frac{{(SR(x,t) – SR_{i,min } )}}{{(SR_{i,max } – SR_{i,min } )}}*(FPAR_{max } – FPAR_{min } ) + FPAR_{min }$$
(15)

where the values of FPARmin and FPARmax are independent of vegetation type and are 0.001 and 0.95, respectively; SRi,max and SRi,min correspond to the 95% and 5% percentiles, respectively, of the ith NDVI. SR(x,t) is represented by the following formula:$$SR(x,t) = frac{1 + NDVI(x,t)}{{1 – NDVI(x,t)}}$$
(16)
A comparison of the estimated results of FPAR-NDVI and FPAR-SR shows that the FPAR estimated by NDVI is higher than the measured value, while the FPAR estimated by SR is lower than the measured value, but the error is less than that estimated directly by NDVI. As a result, these two values can be combined, and their weighted average value is taken as an estimate of the estimated FPAR, while ɑ means weight:$$FPAR(x,t) = alpha FPAR_{NDVI} + (1 – alpha )FPAR_{SR}$$
(17)
Light use efficiency (LUE) refers to the ratio of chemical energy contained in organic dry matter produced per unit area over a certain period of time to the photosynthetically active radiation absorbed by plants projected onto the same area at the same time. Different vegetation types and the same types of vegetation have different light energy utilization rates in different living environments43. The differences are mainly due to the characteristics of the vegetation itself, temperature, moisture, and soil44. Vegetation has the highest utilization rate of light energy under ideal conditions, but the maximum light energy utilization rate in the real environment is mainly affected by temperature and moisture, which can be expressed as follows:$$varepsilon left( {x,t} right) = T_{varepsilon 1} left( {x,t} right) cdot T_{varepsilon 2} left( {x,t} right) cdot W_{varepsilon } left( {x,t} right) cdot varepsilon_{max }$$
(18)
where Tε1(x,t) and Tε2(x,t) represent the stress effects of low temperature and high temperature on light energy utilization, respectively, Wε(x,t) is the effect of water stress on the maximum light energy utilization under ideal conditions, and εmax is the maximum light energy utilization under ideal conditions (gC * MJ−1). The maximum solar energy utilization rate εmax varies depending on the vegetation type. In this study, the maximum light energy utilization rate of different land use types simulated by an improved Carnegie-Ames-Stanford Approach (CASA) model is used as the input parameter of light energy utilization in the CASA model (Table 2). The monthly maximum light energy utilization rate is determined in three steps: first, calculate the APAR, temperature, and water stress factors of all pixels; then, select the NPP measured data of the same time period in the study area; finally, simulate the εmax of vegetation according to the principle of minimum error45. Figure 5 shows the calculation process of NPP. The weight of each habitat type in the HQ is shown in Table 3. The weight value is derived from the official document30. To facilitate the calculation, this paper normalizes the calculation results from 0 to 1 (Fig. 6a).Table 2 Maximum LUE rates of different land use types.Full size tableFigure 5NPP calculation process.Full size imageTable 3 Weight of each habitat type in the HQ.Full size tableVegetation coverage indexThe vegetation coverage index was obtained from the NDVI, which is a simple, effective, and empirical measure of surface vegetation status. The vegetation index mainly describes the difference between the reflection of vegetation in the visible and near-infrared bands and the soil background. This index also reduces the solar elevation angle and noise caused by the atmosphere and is thus the most widely used and effective calculation method. Each vegetation index can be used to quantitatively describe the growth of vegetation under certain conditions. The expression is as follows:$$NDVI = frac{NIR – R}{{NIR + R}}$$
(19)

where NIR and R are reflectance values in the near-infrared and red bands, respectively.NDVI values are obtained by processing the RS images of the Landsat 8 satellite. This satellite is equipped with an operational land imager (OLI) that includes nine bands with a spatial resolution of 30 m, including a 15-m panchromatic band. To facilitate the calculation, this paper normalizes the calculation results from 0 to 1 (Fig. 6b).Figure 6Eco-environment assessment indexes and evaluation rating map (the first quarter was used as an example). (a) Biological abundance index; (b) vegetation coverage index; (c) river density index; (d) land stress index; (e) pollution loading index; and (f) environmental status classification. The Figure is created using ArcGIS ver.10.3 (https://www.esri.com/).Full size imageRiver density indexThe river density index refers to the total length of rivers, lakes, and water resources in the assessed area as a percentage of the assessed area, which is used to reflect the abundance of water in the assessed area and is calculated as follows:$$begin{gathered} River ; density ; index = (A_{riv} *river ; length / area + A_{lak} *water ; area/area hfill \ + A_{res} {*}amount ; of ; resources/area , )/3 hfill \ end{gathered}$$
(20)

where Ariv is the normalization coefficient of river length, with a reference value of 84.3704, Alak is the normalization coefficient of the lake area, with a reference value of 591.7909, and Ares is the normalization coefficient of water resources, with a reference value of 86.387. Finally, the calculation results were normalized from 0 to 1 (Fig. 6c).Land stress indexThe land stress index is the degree to which the land quality in the assessment area is under stress. The weight of the land stress index evaluation is shown in Table 4.Table 4 Weight of the land stress index evaluation.Full size tableThe calculation method is as follows:$$begin{gathered} Land ; stress ; index = A_{ero} *(0.4*severe ; erosion ; area + 0.2*{text{mod}}erate ; erosion ; area hfill \ + 0.2*construction ; land ; area + 0.2*other ; land ; stress ; area)/area hfill \ end{gathered}$$
(21)

where Aero is the normalization coefficient of the land stress index, with a reference value of 236.0436. According to the “Classification criteria for soil erosion”46, the influencing factors of soil erosion, vegetation, soil texture, landform, and precipitation are ranked according to importance. In the calculation of the land stress index, all the land is divided into three categories, in which the weight of severe erosion is 0.4, the weight of non-erosion is 0, and the other erosion types such as moderate erosion and construction land are 0.2. The areas with severe erosion include vegetation coverage less than 30% and areas of soil erosion greater than 3.7 mm/a due to human activities. These areas are generally developed on highly erosive-sensitive soils. Cinnamon soil and loess soil in the study area are highly erosive-sensitive soils. Therefore, the industrial and mining areas of cinnamon and loess soil types are regarded as severe erosion areas. Areas with vegetation coverage greater than 50% are non-eroded areas, so water bodies and woodlands are divided into non-erodible areas. All areas except these two types have a weight of 0.2. Finally, the calculation results were normalized from 0 to 1 (Fig. 6d).Pollution loading indexThe pollution loading index refers to the load of pollutants in a certain area or an environmental element. In this study, the AQI was used to calculate the pollution loading index, and the results were normalized from 0 to 1 (Fig. 6e).The eco-environmental evaluation score was calculated based on the national environmental protection standard according to the weight of each indicator (Fig. 6f).Improved evaluation system and intelligent evaluation modelImproved evaluation systemConsidering that the evaluation factors in the national environmental protection standards are applicable to ordinary areas, areas affected by mines should have more evaluation factors than those in the standards. Thus, an improved evaluation system was proposed. The improved evaluation system has added factors that affect the environment of the mine based on the factors of the original system. The impact of the mine on the environment is reflected in the pollution of the atmosphere, such as dust from open pits and industrial waste from concentrators; the occupation of land by solid waste, such as ore piles and coal piles; soil pollution, such as the diffusion of heavy metals from coal piles, coal mine concentrator plants, and mines; and the increased likelihood of geological disasters, such as collapse caused by underground mining, spontaneous coal combustion and landslides caused by open-pit mining surfaces. Therefore, the improved evaluation system adds an air pollution range, a solid waste area, a geologic hazard range, and a metallic and non-metallic mine soil pollution buffer to the national environmental protection standards.The area of air pollution in mining areas is generally near open pits and concentrator plants. Therefore, the air pollution range was selected within 50 m around the open pits and the concentrator plants. Due to the dilution and dispersion of the air itself, an estimate of the pollution is the reciprocal of the wind speed (Fig. 7a).Figure 7New factors in the improved evaluation system. (a) Air pollution range; (b) solid waste area; (c) geological hazard range; (d1) non-metallic mine soil pollution buffer; and (d2) metal mine soil pollution buffer. The Figure is created using ArcGIS ver.10.3 (https://www.esri.com/).Full size imageMine solid waste pollution includes a large amount of waste rock from open-pit mining and pit mining, coal gangue produced by coal mining, tailings from beneficiation and slag from smelting. These solid wastes are generally piled up near the mining area. They not only occupy large areas of land and induce geological disasters such as landslides and mudslides but also cause chemical pollution, spontaneous combustion, and radiation from radioactive materials due to long-term stacking. This may affect the health and safety of humans and other biological organisms. The scope of the solid waste area is determined by the ore piles, coal piles, and dumping sites (Fig. 7b).Mine geological disasters are caused by a large number of mining wells and rock and soil deformation, as well as serious changes in the geological, hydrogeological, and natural environments of the mining area, endangering human life and property and destroying mining engineering equipment and mining resources. In this study, the geologic hazard range consists of areas with underground mining stopes and coal piles (Fig. 7c).After the pollutants generated by the mining operation enter the soil, physical and mechanical absorption, retention, colloidal physicochemical adsorption, chemical precipitation, bioabsorption, etc. of the suspended pollutants through the soil continue to accumulate in the upper soil. When pollutants reach a certain maximum, they cause deterioration of the soil composition, cycle, properties, and functions and begin to accumulate in plants, which affects the normal growth and development of plants, decreases crop yield and quality, and ultimately affects human health. Metallic and non-metallic minerals have different effects on soil pollution. The pollution of soil by non-metallic minerals mainly occurs in coal mines and coal piles, and the buffer zone is centred on the coal mines and coal piles. Coal production activities can cause heavy metals in coal piles to enter the soil and cause pollution. Due to different types of heavy metals, the range of soil contamination is different47. Combining the non-metallic mineral industrial squares and coal mine-based non-metallic minerals around the heavy metal soil pollution range, the buffers are graded at 30, 200, and 1000 m (Fig. 7d1)43. The metal mines in Gongyi are mainly aluminium ore and iron ore. Referencing the spread range of heavy metal pollution in the soil of aluminium ore and iron ore mines48,49,50,51, the buffers are graded at 50, 100, 300, and 500 m (Fig. 7d2).The four new elements in the improved evaluation system are normalized from 0 to 1 during the calculation.Intelligent evaluation modelArtificial neural networks, decision trees, and SVMs were calculated using IBM SPSS Modeler software to find an intelligent model suitable for environmental assessment of the mine in the study area. Then, several models with high evaluation accuracy were selected. The SVM, CART, and C5.0 models were chosen for further comparison. The sampling points were selected randomly; 700 sampling points were selected from the area away from the mining area; 100 sampling points were selected from the mining area after random sampling by mine type, and these points were used as training samples. Non-mining evaluation scores were based on the national environmental protection standard, while the mining area scores were based on field investigation. In the field investigation, preliminary scoring of the sampling points was conducted according to mining type, mining intensity, air quality, and surrounding environment. Then, a photo of the field was taken at every sampling point in the mine, and experts were invited to further score the area according to the photo. This score is the relative score obtained by referencing the national environmental protection standard.The index layers of the training samples were used as input, and the scores were used as the output to train the machine learning models. The trained models were applied to the entire study area, and all points except the training sample points were used for verification. After further comparison with SVM, CART, and C5.0, the evaluation accuracy rates of the three methods in the mining area and non-mining area were obtained. In the non-mining area, the model evaluation results of various land use types were compared with the national environmental protection standards. The accuracy in various land use types is shown in Table 5. In the mining area, the model evaluation score is compared with the score from the experts, and the obtained accuracy table is shown in Table 6.Table 5 Accuracy of each algorithm in various land use types in non-mining areas.Full size tableTable 6 Accuracy of each algorithm in the mining area.Full size tableIn non-mining areas, the accuracy of the SVM model is significantly better than that of the other two methods. However, in the mining area, the accuracy of the CART model is higher. Therefore, the SVM model was used to evaluate the area away from the mine, and the CART model was used to evaluate the mining area. The evaluation results of these two models were combined to obtain the evaluation map of the entire study area. More

1.Foley, J. A. et al. Global consequences of land use. Science 309, 570–574 (2005).CAS 

Google Scholar 
2.Carpenter, S. R. et al. Science for managing ecosystem services: beyond the Millennium Ecosystem Assessment. Proc. Natl Acad. Sci. USA 106, 1305–1312 (2009).CAS 

Google Scholar 
3.Nelson, E. et al. Modeling multiple ecosystem services, biodiversity conservation, commodity production, and tradeoffs at landscape scales. Front. Ecol. Environ. 7, 4–11 (2009).
Google Scholar 
4.Maes, J. et al. An indicator framework for assessing ecosystem services in support of the EU Biodiversity Strategy to 2020. Ecosyst. Serv. 17, 14–23 (2016).
Google Scholar 
5.Maes, J. et al. Mapping ecosystem services for policy support and decision making in the European Union. Ecosyst. Serv. 1, 31–39 (2012).
Google Scholar 
6.Millennium Ecosystem Assessment. Ecosystems and Human Well-Being: Synthesis (Island Press, 2005).7.Summary for Policymakers. In Global Assessment Report on Biodiversity and Ecosystem Services (IPBES, 2019).8.Martínez-Harms, M. J. & Balvanera, P. Methods for mapping ecosystem service supply: a review. Int. J. Biodivers. Sci. Ecosyst. Serv. Manage. 8, 17–25 (2012).
Google Scholar 
9.Hauck, J. et al. ‘Maps have an air of authority’: potential benefits and challenges of ecosystem service maps at different levels of decision making. Ecosyst. Serv. 4, 25–32 (2013).
Google Scholar 
10.Balvanera, P. et al. Conserving biodiversity and ecosystem services. Science 291, 2047 (2001).CAS 

Google Scholar 
11.Dick, J., Maes, J., Smith, R. I., Paracchini, M. L. & Zulian, G. Cross-scale analysis of ecosystem services identified and assessed at local and European level. Ecol. Indic. 38, 20–30 (2014).
Google Scholar 
12.UK National Ecosystem Assessment. The UK National Ecosystem Assessment Technical Report. (UNEP-WCMC, 2011); http://uknea.unep-wcmc.org/13.Orsi, F., Ciolli, M., Primmer, E., Varumo, L. & Geneletti, D. Mapping hotspots and bundles of forest ecosystem services across the European Union. Land Use Policy 99, 104840 (2020).
Google Scholar 
14.Holland, R. A. et al. The influence of temporal variation on relationships between ecosystem services. Biodivers. Conserv. 20, 3285–3294 (2011).
Google Scholar 
15.Renard, D., Rhemtull, J. M. & Bennett, E. M. Historical dynamics in ecosystem service bundles. Proc. Natl Acad. Sci. USA 112, 13411–13416 (2015).CAS 

Google Scholar 
16.Rukundo, E. et al. Spatio-temporal dynamics of critical ecosystem services in response to agricultural expansion in Rwanda, East Africa. Ecol. Indic. 89, 696–705 (2018).
Google Scholar 
17.Stürck, J., Schulp, C. J. E. & Verburg, P. H. Spatio-temporal dynamics of regulating ecosystem services in Europe—the role of past and future land use change. Appl. Geogr. 63, 121–135 (2015).
Google Scholar 
18.Rau, A. L. et al. Temporal patterns in ecosystem services research: a review and three recommendations. Ambio 49, 1377–1393 (2020).
Google Scholar 
19.Sutherland, I. J., Bennett, E. M. & Gergel, S. E. Recovery trends for multiple ecosystem services reveal non-linear responses and long-term tradeoffs from temperate forest harvesting. For. Ecol. Manage. 374, 61–70 (2016).
Google Scholar 
20.Hansen, M. C., Stehman, S. V. & Potapov, P. V. Quantification of global gross forest cover loss. Proc. Natl Acad. Sci. USA 107, 8650–8655 (2010).CAS 

Google Scholar 
21.Vanhanen, H. et al. Making Boreal Forests Work for People and Nature (IUFRO, 2012).22.Pan, Y. et al. A large and persistent carbon sink in the world’s forests. Science 333, 988–993 (2011).CAS 

Google Scholar 
23.Moen, J. et al. Eye on the Taiga: removing global policy impediments to safeguard the boreal forest. Conserv. Lett. 7, 408–418 (2014).
Google Scholar 
24.Global Forest Industry (Swedish Forest Industries, 2019); https://www.forestindustries.se/forest-industry/statistics/global-forest-industry/25.Saastamoinen, O., Kangas, K. & Aho, H. The picking of wild berries in Finland in 1997 and 1998. Scand. J. For. Res. 15, 645–650 (2000).
Google Scholar 
26.Gamfeldt, L. et al. Higher levels of multiple ecosystem services are found in forests with more tree species. Nat. Commun. 4, 1340 (2013).
Google Scholar 
27.Hou, Y., Li, B., Müller, F., Fu, Q. & Chen, W. A conservation decision-making framework based on ecosystem service hotspot and interaction analyses on multiple scales. Sci. Total Environ. 643, 277–291 (2018).CAS 

Google Scholar 
28.Blumstein, M. & Thompson, J. R. Land-use impacts on the quantity and configuration of ecosystem service provisioning in Massachusetts, USA. J. Appl. Ecol. 52, 1009–1019 (2015).
Google Scholar 
29.Fernandez-Campo, M., Rodríguez-Morales, B., Dramstad, W. E., Fjellstad, W. & Diaz-Varela, E. R. Ecosystem services mapping for detection of bundles, synergies and trade-offs: examples from two Norwegian municipalities. Ecosyst. Serv. 28, 283–297 (2017).
Google Scholar 
30.Gissi, E., Fraschetti, S. & Micheli, F. Incorporating change in marine spatial planning: a review. Environ. Sci. Policy 92, 191–200 (2019).
Google Scholar 
31.Maxwell, S. M., Gjerde, K. M., Conners, M. G. & Crowder, L. B. Mobile protected areas for biodiversity on the high seas. Science 367, 252–254 (2020).CAS 

Google Scholar 
32.Willcock, S. et al. Do ecosystem service maps and models meet stakeholders’ needs? A preliminary survey across sub-Saharan Africa. Ecosyst. Serv. 18, 110–117 (2016).
Google Scholar 
33.Jonsson, M., Bengtsson, J., Gamfeldt, L., Moen, J. & Snäll, T. Levels of forest ecosystem services depend on specific mixtures of commercial tree species. Nat. Plants 5, 141–147 (2019).
Google Scholar 
34.Pohjanmies, T. et al. Impacts of forestry on boreal forests: an ecosystem services perspective. Ambio 46, 743–755 (2017).
Google Scholar 
35.Miina, J., Hotanen, J.-P. & Salo, K. Modelling the abundance and temporal variation in the production of bilberry (Vaccinium myrtillus L.) in Finnish mineral soil forests. Silva Fenn. 43, 577–593 (2009).
Google Scholar 
36.Hertel, A. G. et al. Berry production drives bottom–up effects on body mass and reproductive success in an omnivore. Oikos 127, 197–207 (2018).
Google Scholar 
37.Thiffault, E. Boreal forests and soils. Dev. Soil Sci. 36, 59–82 (2019).
Google Scholar 
38.Jonsson, M., Bengtsson, J., Moen, J., Gamfeldt, L. & Snäll, T. Stand age and climate influence forest ecosystem service delivery and multifunctionality. Environ. Res. Lett. 15, 0940a8 (2020).
Google Scholar 
39.Stokland, J. N. Volume increment and carbon dynamics in boreal forest when extending the rotation length towards biologically old stands. For. Ecol. Manage. 488, 119017 (2021).
Google Scholar 
40.Harmon, M. E., Ferrell, W. K. & Franklin, J. F. Effects on carbon storage of conversion of old-growth forests to young forests. Science 247, 699–702 (1990).CAS 

Google Scholar 
41.Mazziotta, A. et al. Applying a framework for landscape planning under climate change for the conservation of biodiversity in the Finnish boreal forest. Glob. Change Biol. 21, 637–651 (2015).
Google Scholar 
42.Triviño, M. et al. Optimizing management to enhance multifunctionality in a boreal forest landscape. J. Appl. Ecol. 54, 61–70 (2017).
Google Scholar 
43.Qiu, J. & Turner, M. G. Spatial interactions among ecosystem services in an urbanizing agricultural watershed. Proc. Natl Acad. Sci. USA 110, 12149–12154 (2013).CAS 

Google Scholar 
44.Felipe-Lucia, M. R. et al. Multiple forest attributes underpin the supply of multiple ecosystem services. Nat. Commun. 9, 4839 (2018).
Google Scholar 
45.Eggers, J., Räty, M., Öhman, K. & Snäll, T. How well do stakeholder-defined forest management scenarios balance economic and ecological forest values? Forests 11, 86 (2020).
Google Scholar 
46.Eyvindson, K., Repo, A. & Mönkkönen, M. Mitigating forest biodiversity and ecosystem service losses in the era of bio-based economy. For. Policy Econ. 92, 119–127 (2018).
Google Scholar 
47.Rusch, A., Bommarco, R., Jonsson, M., Smith, H. G. & Ekbom, B. Flow and stability of natural pest control services depend on complexity and crop rotation at the landscape scale. J. Appl. Ecol. 50, 345–354 (2013).
Google Scholar 
48.Schipanski, M. E. et al. A framework for evaluating ecosystem services provided by cover crops in agroecosystems. Agric. Syst. 125, 12–22 (2014).
Google Scholar 
49.Hufnagel, J., Reckling, M. & Ewert, F. Diverse approaches to crop diversification in agricultural research. A review. Agron. Sustain. Dev. 40, 14 (2020).
Google Scholar 
50.Guerry, A. D. et al. Modeling benefits from nature: using ecosystem services to inform coastal and marine spatial planning. Int. J. Biodivers. Sci. Ecosyst. Serv. Manage. 8, 107–121 (2012).
Google Scholar 
51.Wikström, P. et al. The Heureka Forestry Decision Support System: An Overview. Math. Comput. For Nat.-Resour. Sci. 3, 87–94 (2011).
Google Scholar 
52.Forest Statistics (Swedish University of Agricultural Sciences, 2020).53.Eriksson, A., Snäll, T. & Harrison, P. J. Analys av miljöförhållanden ‐ SKA 15. Report 11 (Swedish Forest Agency, 2015).54.Axelsson, A.-L. et al. in National Forest Inventories—Pathways for Common Reporting (eds Tomppo, E. et al.) 541–553 (Springer, 2010).55.Marklund, L. G. Biomass Functions for Pine, Spruce and Birch in Sweden (1988).56.Petersson, H. & Ståhl, G. Functions for below-ground biomass of Pinus sylvestris, Picea abies, Betula pendula and Betula pubescens in Sweden. Scand. J. For. Res. 21, 24–83 (2006).
Google Scholar 
57.Miina, J., Pukkala, T. & Kurttila, M. Optimal multi-product management of stands producing timber and wild berries. Eur. J. For. Res. 135, 781–794 (2016).
Google Scholar 
58.Schröter, M. & Remme, R. P. Spatial prioritisation for conserving ecosystem services: comparing hotspots with heuristic optimisation. Landsc. Ecol. 31, 431–450 (2016).
Google Scholar 
59.Wu, J., Feng, Z., Gao, Y. & Peng, J. Hotspot and relationship identification in multiple landscape services: a case study on an area with intensive human activities. Ecol. Indic. 29, 529–537 (2013).CAS 

Google Scholar 
60.Akaike, H. A new look at the statistical model identification. IEEE Trans. Automat. Contr. 19, 716–723 (1974).
Google Scholar 
61.R: A Language and Environment for Statistical Computing (R Development Core Team, 2014); https://www.R-project.org/ More