Philippa Kaur

1.Arrhenius, O. Species and area. J. Ecol. 9, 95–99 (1921).
Google Scholar 
2.MacArthur, R. H. & Wilson, E. O. The Theory of Island Biogeography (Princeton Univ. Press, 1967).3.Rosenzweig, M. L. Species Diversity in Space and Time (Cambridge Univ. Press, 1995).4.Smith, A. B., Sandel, B., Kraft, N. J. B. & Carey, S. Characterizing scale‐dependent community assembly using the functional‐diversity–area relationship. Ecology 94, 2392–2402 (2013).PubMed 

Google Scholar 
5.Mazel, F. et al. Multifaceted diversity–area relationships reveal global hotspots of mammalian species, trait and lineage diversity. Glob. Ecol. Biogeogr. 23, 836–847 (2014).PubMed 
PubMed Central 

Google Scholar 
6.Dias, R. A. et al. Species richness and patterns of overdispersion, clustering and randomness shape phylogenetic and functional diversity–area relationships in habitat islands. J. Biogeogr. 47, 1638–1648 (2020).
Google Scholar 
7.Pereira, H. M. et al. Scenarios for global biodiversity in the 21st century. Science 330, 1496–1501 (2010).CAS 
PubMed 

Google Scholar 
8.Pimm, S. L., Russell, G. J., Gittleman, J. L. & Brooks, T. M. The future of biodiversity. Science 269, 347–350 (1995).CAS 
PubMed 

Google Scholar 
9.Simberloff, D. in Tropical Deforestation and Species Extinction (eds Whitmore, T. C. & Sayer, J. A.) 75–89 (Chapman & Hall, 1992).10.Jordano, P. Chasing ecological interactions. PLoS Biol. 14, e1002559 (2016).PubMed 
PubMed Central 

Google Scholar 
11.Montoya, J. M., Woodward, G., Emmerson, M. C. & Solé, R. V. Press perturbations and indirect effects in real food webs. Ecology 90, 2426–2433 (2009).PubMed 

Google Scholar 
12.Lurgi, M., López, B. C., Montoya, J. M. & Lopez, B. C. Novel communities from climate change. Philos. Trans. R. Soc. Lond. B 367, 2913–2922 (2012).
Google Scholar 
13.Tylianakis, J. M., Tscharntke, T. & Lewis, O. T. Habitat modification alters the structure of tropical host–parasitoid food webs. Nature 445, 202–205 (2007).CAS 
PubMed 

Google Scholar 
14.Montoya, J. M., Rodriguez, M. Á. & Hawkins, B. A. Food web complexity and higher-level ecosystem services. Ecol. Lett. 6, 587–593 (2003).
Google Scholar 
15.Reiss, J., Bridle, J. R., Montoya, J. M. & Woodward, G. Emerging horizons in biodiversity and ecosystem functioning research. Trends Ecol. Evol. 24, 505–514 (2009).PubMed 

Google Scholar 
16.Thompson, R. M. et al. Food webs: reconciling the structure and function of biodiversity. Trends Ecol. Evol. 27, 689–697 (2012).PubMed 

Google Scholar 
17.Cohen, J. E. & Newman, C. M. Community area and food-chain length: theoretical predictions. Am. Nat. 138, 1542–1554 (1991).
Google Scholar 
18.Schoener, T. W. Food webs from the small to the large: the Robert H. MacArthur Award lecture. Ecology 70, 1559–1589 (1989).
Google Scholar 
19.Post, D. M., Pace, M. L. & Hairston, N. G. Ecosystem size determines food-chain length in lakes. Nature 405, 1047–1049 (2000).CAS 
PubMed 

Google Scholar 
20.Brose, U., Ostling, A., Harrison, K. & Martinez, N. D. Unified spatial scaling of species and their trophic interactions. Nature 428, 167–171 (2004).CAS 
PubMed 

Google Scholar 
21.Galiana, N. et al. The spatial scaling of species interaction networks. Nat. Ecol. Evol. 2, 782–790 (2018).PubMed 

Google Scholar 
22.Wood, S. A., Russell, R., Hanson, D., Williams, R. J. & Dunne, J. A. Effects of spatial scale of sampling on food web structure. Ecol. Evol. 5, 3769–3782 (2015).PubMed 
PubMed Central 

Google Scholar 
23.Pimm, S. L. et al. Food web patterns and their consequences. Nature 350, 669–674 (1991).
Google Scholar 
24.Martinez, N. D. Constant connectance in community food webs. Am. Nat. 139, 1208–1218 (1992).
Google Scholar 
25.Ings, T. C. et al. Ecological networks–beyond food webs. J. Anim. Ecol. 78, 253–69 (2009).PubMed 

Google Scholar 
26.Montoya, J. M. & Solé, R. V. Topological properties of food webs: from real data to community assembly models. Oikos 102, 614–622 (2003).
Google Scholar 
27.Drakare, S., Lennon, J. J. & Hillebrand, H. The imprint of the geographical, evolutionary and ecological context on species–area relationships. Ecol. Lett. 9, 215–227 (2006).PubMed 

Google Scholar 
28.Preston, F. W. Time and space and the variation of species. Ecology 41, 611–627 (1960).
Google Scholar 
29.Turner, W. R. & Tjørve, E. Scale-dependence in species–area relationships. Ecography 6, 721–730 (2005).
Google Scholar 
30.Bengtsson, J. Confounding variables and independent observations in comparative analyses of food webs. Ecology 75, 1282–1288 (1994).
Google Scholar 
31.Vermaat, J. E., Dunne, J. A. & Gilbert, A. J. Major dimensions in food-web structure properties. Ecology 90, 278–282 (2009).PubMed 

Google Scholar 
32.Dunne, J. A. et al. Parasites affect food web structure primarily through increased diversity and complexity. PLoS Biol. 11, e1001579 (2013).CAS 
PubMed 
PubMed Central 

Google Scholar 
33.Poisot, T. & Gravel, D. When is an ecological network complex? Connectance drives degree distribution and emerging network properties. PeerJ 2, e251 (2014).PubMed 
PubMed Central 

Google Scholar 
34.Cohen, J. E. & Briand, Fredeiri Trophic links of community food webs. Proc. Natl Acad. Sci. USA 81, 4105–4109 (1984).CAS 
PubMed 
PubMed Central 

Google Scholar 
35.Roslin, T., Várkonyi, G., Koponen, M., Vikberg, V. & Nieminen, M. Species–area relationships across four trophic levels—decreasing island size truncates food chains. Ecography 37, 443–453 (2014).
Google Scholar 
36.Holt, R. D., Lawton, J. H., Polis, G. A. & Martinez, N. D. Trophic rank and the species–area relationship. Ecology 80, 1495–1504 (1999).
Google Scholar 
37.Dunne, J. A., Williams, R. J. & Martinez, N. D. Food-web structure and network theory: the role of connectance and size. Proc. Natl Acad. Sci. USA 99, 12917–12922 (2002).CAS 
PubMed 
PubMed Central 

Google Scholar 
38.Montoya, J. M., Pimm, S. L. & Solé, R. V. Ecological networks and their fragility. Nature 442, 259–264 (2006).CAS 
PubMed 

Google Scholar 
39.James, A., Pitchford, J. W. & Plank, M. J. Disentangling nestedness from models of ecological complexity. Nature 487, 227–230 (2012).CAS 
PubMed 

Google Scholar 
40.Valverde, S. et al. The architecture of mutualistic networks as an evolutionary spandrel. Nat. Ecol. Evol. 2, 94–99 (2018).PubMed 

Google Scholar 
41.Valiente-Banuet, A. et al. Beyond species loss: the extinction of ecological interactions in a changing world. Funct. Ecol. 29, 299–307 (2015).
Google Scholar 
42.Janzen, D. H. The deflowering of central America. Nat. Hist. 83, 49–53 (1974).43.Mendoza, M. & Araújo, M. B. Climate shapes mammal community trophic structures and humans simplify them. Nat. Commun. 10, 5197 (2019).PubMed 
PubMed Central 

Google Scholar 
44.Emer, C. et al. Seed dispersal networks in tropical forest fragments: area effects, remnant species, and interaction diversity. Biotropica 52, 81–89 (2020).
Google Scholar 
45.McWilliams, C., Lurgi, M., Montoya, J. M., Sauve, A. & Montoya, D. The stability of multitrophic communities under habitat loss. Nat. Commun. 10, 2322 (2019).PubMed 
PubMed Central 

Google Scholar 
46.McCann, K. S. The diversity–stability debate. Nature 405, 228–233 (2000).CAS 
PubMed 

Google Scholar 
47.Fig, T., Mccann, K., Hastings, A. & Huxel, G. R. Weak trophic interactions and the balance of nature. Nature 395, 794–798 (1998).
Google Scholar 
48.Pimm, S. L. & Lawton, J. H. Are food webs divided into compartments? J. Anim. Ecol. 49, 879–898 (1980).
Google Scholar 
49.Macfadyen, S., Gibson, R. H., Symondson, W. O. C. & Memmott, J. Landscape structure influences modularity patterns in farm food webs: consequences for pest control. Ecol. Appl. 21, 516–524 (2011).PubMed 
PubMed Central 

Google Scholar 
50.Reverté, S. et al. Spatial variability in a plant–pollinator community across a continuous habitat: high heterogeneity in the face of apparent uniformity. Ecography 42, 1558–1568 (2019).
Google Scholar 
51.Torné‐Noguera, A., Arnan, X., Rodrigo, A. & Bosch, J. Spatial variability of hosts, parasitoids and their interactions across a homogeneous landscape. Ecol. Evol. 10, 3696–3705 (2020).PubMed 
PubMed Central 

Google Scholar 
52.Hernández‐Castellano, C. et al. A new native plant in the neighborhood: effects on plant–pollinator networks, pollination, and plant reproductive success. Ecology 101, e03046 (2020).PubMed 

Google Scholar 
53.Osorio, S., Arnan, X., Bassols, E., Vicens, N. & Bosch, J. Local and landscape effects in a host–parasitoid interaction network along a forest–cropland gradient. Ecol. Appl. 25, 1869–1879 (2015).PubMed 

Google Scholar 
54.Kaartinen, R. & Roslin, T. Shrinking by numbers: landscape context affects the species composition but not the quantitative structure of local food webs. J. Anim. Ecol. 80, 622–631 (2011).PubMed 

Google Scholar 
55.Vázquez, D. P. & Simberloff, D. Changes in interaction biodiversity induced by an introduced ungulate. Ecol. Lett. 6, 1077–1083 (2003).
Google Scholar 
56.Mulder, C., Den Hollander, H. A. & Hendriks, A. J. Aboveground herbivory shapes the biomass distribution and flux of soil invertebrates. PLoS ONE 3, e3573 (2008).PubMed 
PubMed Central 

Google Scholar 
57.Montoya, D., Yallop, M. L. & Memmott, J. Functional group diversity increases with modularity in complex food webs. Nat. Commun. 6, 7379 (2015).CAS 
PubMed 

Google Scholar 
58.Grass, I., Jauker, B., Steffan-Dewenter, I., Tscharntke, T. & Jauker, F. Past and potential future effects of habitat fragmentation on structure and stability of plant–pollinator and host–parasitoid networks. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-018-0631-2 (2018).59.Cagnolo, L., Salvo, A. & Valladares, G. Network topology: patterns and mechanisms in plant–herbivore and host–parasitoid food webs. J. Anim. Ecol. 80, 342–351 (2011).PubMed 

Google Scholar 
60.Maiorano, L., Montemaggiori, A., Ficetola, G. F., O’Connor, L. & Thuiller, W. TETRA‐EU 1.0: a species‐level trophic metaweb of European tetrapods. Glob. Ecol. Biogeogr. 29, 1452–1457 (2020).61.Kopelke, J. et al. Food‐web structure of willow‐galling sawflies and their natural enemies across Europe. Ecology 98, 1730 (2017).PubMed 

Google Scholar 
62.Sole, R. V. & Montoya, M. Complexity and fragility in ecological networks. Proc. R. Soc. Lond. B 268, 2039–2045 (2001).CAS 

Google Scholar 
63.Broido, A. D. & Clauset, A. Scale-free networks are rare. Nat. Commun. 10, 1017 (2019).PubMed 
PubMed Central 

Google Scholar 
64.Guilhaumon, F., Mouillot, D. & Gimenez, O. mmSAR: an R-package for multimodel species–area relationship inference. Ecography 33, 420–424 (2010).
Google Scholar 
65.Matthews, T. J., Triantis, K. A., Whittaker, R. J. & Guilhaumon, F. sars: an R package for fitting, evaluating and comparing species–area relationship models. Ecography https://doi.org/10.1111/ecog.04271 (2019).66.Galiana, N. Ecological network complexity scales with area. Dryad https://doi.org/10.5061/dryad.zcrjdfndg (2021). More

Hydrothermal vent rock and sediment sample collectionRock no. NA091-R045 (source of Ca. H. endolithica PR6, Ca. H. repetitus FW102 and Thorarchaeote FW25) and rock no. NA091-R008 (source of Heimdall group Gerdarchaeote AC18) were retrieved from the Auka hydrothermal vent site situated on the margin of the southern Pescadero Basin of the Gulf of California using remotely operated vehicle Hercules during research expedition NA091 on E/V Nautilus on 2 November 2017. Local venting fluids have a measured temperature approaching 300 °C, contain hydrocarbons and hydrogen and are precipitating minerals, such as calcite and barite15. R045 was collected during dive H1658 at coordinates 23.956987786° N, 108.86227922° W at a water depth of 3,674 m, near shimmering water, a sign of locally focused hydrothermal fluid discharge. R008 was collected during dive H1657 at coordinates 23° 57′ N, 108° 52′ W at a water depth of 3,651 m. After shipboard recovery, rock samples were placed in Mylar bags prefilled with 0.2 µm filtered bottom seawater collected during the same dive, flushed with N2 gas for 10 min, sealed and stored at 4 °C until preparation for incubations in the laboratory.Sediment sample no. FK181031-S0193-PC3 (source of Ca. H. aukensis) was collected during the research expedition FK181031 on R/V Falkor to the southern Pescadero Basin on 14 November 2018. The sample was collected during dive S193 at the Auka hydrothermal vent site (23.954822° N, 108.863009° W, water depth of 3,657 m), near the site where rocks nos. NA091-R045 and NA091-R008 were collected in 2017. The sediment push core was extruded upwards and sectioned into discrete 3 cm depth horizons on board immediately after recovery, transferred into sterile Whirl-Pak bags and sealed in a larger Mylar bag, flushed with argon gas, heat-sealed and stored at 4 °C until use in the laboratory.Sample collection permits for the expedition were granted by the Dirección General de Ordenamiento Pesquero y Acuícola, Comisión Nacional de Acuacultura y Pesca (Permiso de Pesca de Fomento no. PPFE/DGOPA-200/18) and the Dirección General de Geografía y Medio Ambiente, Instituto Nacional de Estadística y Geografía (authorization no. EG0122018), with the associated diplomatic note no. 18-2083 (CTC/07345/18) from the Secretaría de Relaciones Exteriores-Agencia Mexicana de Cooperación Internacional para el Desarrollo/Dirección General de Cooperación Técnica y Científica.Artificial seawater medium recipeArtificial seawater was prepared as described in Scheller et al.47 with minor modifications. Briefly, 1 l of artificial seawater (ASW) medium contained 46.6 mM MgCl2, 9.2 mM CaCl2, 485 mM NaCl, 7 mM KCl, 20 mM Na2SO4, 1 mM K2HPO4, 2 mM NH4Cl, 1 ml of 1,000× trace element solution, 1 ml of 1,000× vitamin solution and 0.5 mg of resazurin and was buffered by 25 mM HEPES buffer adjusted to pH 7.5. One litre of 1,000× trace element solution contained 50 mM nitrilotriacetic acid, 5 mM FeCl3, 2.5 mM MnCl2, 1.3 mM CoCl2, 1.5 mM ZnCl2, 0.32 mM H3BO3, 0.38 mM NiCl2, 0.03 mM Na2SeO3, 0.01 mM CuCl2, 0.21 mM Na2MoO4 and 0.02 mM Na2WO4. One litre of 1,000× vitamin solution contained 82 μM d-biotin, 45 μM folic acid, 490 μM pyridoxine, 150 μM thiamine, 410 μM nicotinic acid, 210 μM pantothenic acid, 310 μM para-aminobenzoic acid, 240 μM lipoic acid, 14 μM choline chloride and 7.4 μM vitamin B12.Enrichment cultivationRock no. NA091-R045 was anaerobically fragmented; then, approximately 5 g wet weight was crushed using a sterile agate mortar and pestle on 8 November 2018 and immediately immersed in anaerobic ASW medium in 25–125 ml of butyl rubber-stoppered serum bottles supplemented with different carbon/energy sources, including lactate, H2/CO2, hexane and decane and incubated in the dark at 40 °C (Extended Data Fig. 1a). The headspace for all cultures was flushed and overpressurized with N2 gas (2 atm). For the H2-containing cultures, the N2 gas headspace was replaced with H2/CO2 at an 80:20 mixture by flushing for 1 min and subsequent equilibration at 2 atm. After 33 d of incubation, the lactate-fed first-generation culture produced 5 mM sulphide, indicating active sulphate reduction. This enrichment was mixed by gentle shaking and diluted 1:100 vol/vol into fresh anaerobic ASW medium containing the same suite of carbon/energy sources as described above (Extended Data Fig. 1b). A transfer using the liquid fraction-lacking rock particles from the primary lactate enrichment was also included to enrich for members of the planktonic community alone with lactate as the carbon and energy source. This enrichment was later found to be devoid of the AAG (Heimdall) phylotype. Third- and fourth-generation cultures were set up in the following months through 1:100 dilution (Extended Data Fig. 1b). Further details of microbial community development in these enrichments are provided in Supplementary Note 1 and Supplementary Tables 1–3.R008 was prepared as above except using 2 atm of methane in the headspace as the sole carbon source and electron donor. The culture was passaged twice using a 1:100 dilution under the same culturing conditions; the cell fraction was collected by centrifugation after a total of 22 months for metagenomic sequencing (described below).For sediment enrichment cultivation, the top 3 cm section of the sediment core was mixed with anaerobic ASW at a 1:4 vol/vol ratio; a total of 60 ml volume each was dispensed into seven 125 ml glass serum bottles sealed with butyl rubber stoppers. The headspace was replaced by ethane (2 atm) in 2 bottles (Supplementary Table 5), while the headspace in 1 bottle was replaced by 100% N2 gas (2 atm). The cultures were incubated at 37 °C in the dark. Further details on microbial community development are provided in Supplementary Note 1 and Supplementary Table 4.Mineralogical analysesThe mineralogical composition of rocks NA091-R045 and R008 was characterized on a PANalytical X’Pert Pro X-Ray diffractometer. A dried rock aliquot was finely powdered using a clean agate mortar and pestle and scanned from 3 to 75° (2θ angle) at a 0.0167° step size. Mineral identification was performed with the X’Pert HighScore software v4.1 using the search and march algorithm.DNA extractionCombined cells with rock or sediment substrate were pelleted through centrifugation at 13,000 r.p.m. for 3 min. For amplicon sequencing, unless specified in Supplementary Table 6, DNA was extracted using the Qiagen DNeasy PowerSoil kit (catalogue no. 47014) according to the manufacturer’s instructions as described previously48 with a minor modification, where mechanical shearing was carried out using the MP Biomedicals FastPrep-24 system (catalogue no. 116004500) at level 5.5 for 45 s. For genomic sequencing, incubated rock and sediment cultures were extracted using multiple approaches, including the Qiagen DNeasy PowerSoil kit, ZymoBIOMICS 96 MagBead DNA Kit (catalogue no. D4302; Zymo Research Corporation), Quick-DNA 96 Kit (catalogue no. D3010; Zymo Research Corporation), ZymoBIOMICS DNA Microprep Kit (catalogue no. D4301; Zymo Research Corporation) and a standard phenol/chloroform-based protocol. The list of samples and their extraction methods are provided in Supplementary Table 6.16S rRNA gene amplicon sequencingFor amplicon (iTAG) sequencing of 16S rRNA genes, extracted DNA was amplified using primer pair 515f/806r GTGCCAGCMGCCGCGGTAA/ GGACTACHVGGGTWTCTAAT, barcoded and sequenced at Laragen using the Illumina MiSeq platform and analysed using Qiime v.1.8.0 (ref. 49) as described previously48. Taxonomic assignment was based on the SILVA 138 database (https://www.arb-silva.de)50.Full-length 16S archaeal rRNA gene sequences were amplified using the archaeal primer pair SSU1Arf/SSU1492Rngs TCCGGTTGATCCYGCBRG/ CGGNTACCTTGTKACGAC as described by Bahram et al.51, multiplexed as instructed by PacBio and sequenced using the PacBio Sequel II at the Brigham Young University DNA Sequencing Center and then analysed using the DADA2 package v1.9.1 in R v3.6.0 as described in Callahan et al.52 using the SILVA 138 database for taxonomic classification. Note that in the SILVA 138 database, all Asgard archaea clades are classified under Asgardarchaeota.Metagenomic sequencingA total of 11 metagenomic sequencing runs were performed using the Illumina and Oxford Nanopore platforms, with details listed in Supplementary Table 6. For Illumina short-read sequencing, libraries were constructed using the NEBNext Ultra and Nextera Flex Library kits as specified in the Supplementary Table 6. Sequencing was carried out using a HiSeq 2500 system (single-end, 100 bp) at the Caltech Genetics and Genomics Laboratory and HiSeq 4000 system at Novogene (paired-end, 150 bp). Only paired-end data were used for assembly, while all data were used for error correction. Due to the low DNA quantity obtained from the sediment incubation that yielded Ca. H. aukensis, we used multiple displacement amplification with the QIAGEN REPLI g Midi Kit before library preparation for Nanopore sequencing. Oxford Nanopore sequencing libraries were constructed using the PCR Barcoding Kit (catalogue no. SQK-PBK004) and were sequenced on MinION flow cells FLO-MIN106. Base calling was performed with the ONT Guppy software v.3.4.5.Genome assembly, error correction and read coverage mappingTwo different approaches were used to assemble contiguous genomes from metagenomes. For species of interest, if Nanopore sequencing yielded high read coverage and read lengths N50  > 2 kb, we obtained more contiguous genomes through de novo assembly purely based on Nanopore reads. If Nanopore sequencing did not yield a high number of reads or exhibited low read lengths, we obtained more contiguous genomes through de novo assembly first based on Illumina reads and then joined using Nanopore reads.For Ca. H. endolithica, Nanopore sequencing data were assembled de novo using Canu17 v.2.1, which yielded a 30 Mbp assembly, including a 3.4 Mbp contig. The approximate 40 kilobase (kb) regions at two ends of an approximate 3.4 Mbp contig were repetitive. This repeated region was deleted at one end and the two ends were joined to result in a circular genome. The resulting genome was mapped using BamM (http://ecogenomics.github.io/BamM/, based on Burrows–Wheeler Aligner53 mapping) with 150 bp Illumina paired-end reads (88× coverage on average) and 100 bp single-end reads (20× coverage). Mapped reads were then used for error correction through pilon54 v.1.22. To account for the reduced mapping at the edges (approximate 50 bp region), the two ends of the genomic sequence were joined, read-mapped and error-corrected again using the same methods. After the genome was annotated, it was rotated such that the genomic sequence ended with tRNA (GlyCCC), which was the integration site of the putative provirus HeimV1. All sequencing reads derived from incubations of the same rock were mapped onto the final genome using BamM, which was then used for coverage calculation through bedtools (https://bedtools.readthedocs.io/en/latest/).For Ca. H. aukensis, Illumina PE150 bp sequencing data were assembled using SPAdes18 v.3.14.1 with the ‘-meta’ option and k-mers 21,33,55,77,99. The assembly was then scaffolded using Nanopore reads through two iterations of LRScaf55 v.1.1.10. The Ca. H. aukensis genome was joined after trimming the identical sequences at the two ends. The end-joining region was verified through PCR amplification and Sanger sequencing using the primer pair CGCTTTCTTCAAACAATATTTCTGGTG/CTTACTTTCTCTCGGTCCATTTTTCAC. Finally, a 1 kbp stretch of unresolved genomic sequence at an approximate 2.9 Mbp position was resequenced through PCR amplification and Sanger sequencing using the primers GAGTTTTTTCAATCTTATAATGCCAAACTAAAAAATAG (forward), CAGTCAGATTTGACACAATTTTGGTC (reverse) and GCTGGACTCAACCTATAACTAATAGT (reverse). The final assembly was read-mapped, error-corrected through pilon v.1.24 using 346× coverage. It was rotated as described above to place the tRNA gene GlyCCC at the end.The metagenome containing the Lokiarchaeote Ca. H. repetitus FW102 was assembled using Canu v.2.1, as described for the Ca. H. endolithica genome, and then binned using metabat2 v.2.15 (ref. 56) with default parameters. The bin was then used to recruit long reads using minimap2 v.2.17 and reassembled and binned again. We then used LRScaf to scaffold the contigs and used ten iterations of pilon v.1.24 to achieve error correction and resolve ambiguous bases.The Thorarcheote FW25 MAG was assembled using the hybrid assembly of Illumina reads and Nanopore reads using SPAdes v.3.14.1 with k-mers 21,33,55,77,99, and then binned using metabat2 v.2.15 with default parameters. The MAG bin was then used to recruit reads through MIRAbait in the MIRA v.4 package (http://mira-assembler.sourceforge.net/docs/DefinitiveGuideToMIRA.html#chap_intro). These reads were then used for hybrid assembly with Nanopore long reads via SPAdes v.3.14.1 with k-mers 21,33,55,77,99. It was then binned again using metabat2 v.2.15 with default parameters to yield the final Thorarcheote FW25 MAG.The metagenome containing Gerdarchaeote AC18 was assembled from Illumina reads using SPAdes v.3.14.1 with k-mers 21,33,55,77,99 and then binned using metabat2 v.2.15 with default parameters. The MAG bin was then used to recruit reads through MIRAbait in the MIRA v.4 package and then reassembled and binned using SPAdes and metabat2 to yield the final Gerdarchaeote AC18 bin.Alignment fraction, ANI and AAIANI and alignment fraction values, independently calculated for rRNA, tRNA and coding gene sequences were obtained using ANIcalculator57 2014-127, v.1.0 (https://ani.jgi.doe.gov/html/download.php?). Note that Lokiarchaeote FW102 contains 2 copies of 16S rRNA genes at 99% identity with each other, and Thorarchaeote BC has a partial 16S rRNA gene. The alignment of 16S rRNA was carried out using SINA58 v.1.2.11. The AAI values of translated proteomes were obtained with the enveomics package v1.8.059. The final output is shown in Supplementary Table 7.Genome and mobilome annotationsGene calling was done using a combination of Prodigal v.2.6.3 and Glimmer v.3.0.2 using translation code 11 within the RASTtk60 pipeline, now under the PATRIC package v1.03261. Translated coding sequences were annotated and domain-assigned using eggNOG mapper39 v.2. The tRNA, 16S rRNA and 23S rRNA genes were identified using RNAmmer62 v.1.2 embedded in RASTtk. Thus far, 5S rRNA gene sequences could not be predicted through the existing HMM using various approaches. Long, non-tandem repeats were identified using RASTtk with the default cut-off of 95% identity and 100 bp. Tandem repeat sequences were identified using RASTtk, Prokka v1.14.6 and CRISPRCasTyper 1.1.463. Prokka and CRISPRCasTyper both employ MinCED (https://github.com/ctSkennerton/minced) to identify repeats and detect intragenic tandem repeats, which were manually removed from the CRISPR–Cas analyses. The Cas genes were annotated using CRISRCasTyper.All identified Heimdallarchaeum mobilomes were further analysed using PSI-BLAST 1.10.064, CDD search v3.1965 and PhANNs webserver (version March 2021)37.Genome evaluation and HMM constructionMarker coverage was carried out using a two-step process. First, we used the automated marker analyses via CheckM66 v.1.1.3 with the lineage_wf option and the default HMM E value cut-off, which included the 149 standard archaeal single-copy marker set. Next, each of the missing markers was examined with hmmer67 v.3.3.2 using the hmmsearch option with manual inspection of alignment regions and bitscores. This rescued markers unidentified through the default cut-offs by CheckM as well as divergent variants that most likely functionally replace the genuinely missing marker. The detailed description of markers missed by CheckM can be found in Supplementary Note 2 and the final evaluation of marker presence is displayed in Extended Data Fig. 4a and Supplementary Table 15. Next, we constructed an updated HMM set to replace the CheckM set by (1) updating all HMM to the most recent versions, (2) removing the six commonly missing or duplicated markers shown in Extended Data Fig. 4a from the list and (3) overcoming the pitfall of existing HMMs constructed using only a few sequences acquired from Euryarchaeota and Crenarchaeota. We manually constructed Asgard-specific versions based on the 282 Asgard archaea genomes. The HMMs constructed in this study are PF00832.ASG, PF00861.ASG, PF01194.ASG, PF01287.ASG, PF01667.ASG, PF03874.ASG, PF03876.ASG, PF13656.ASG, TIGR00270.ASG, TIGR00336.ASG, TIGR00442.ASG, TIGR02338.ASG and TIGR03677.ASG. The updated HMM file has been provided as a supplementary data file. The updated HMM was used to evaluate the 282 genomes reported in this study and in the literature3,6,7,8,9,10,11,12,16,23,26,68,69,70,71,72,73,74,75,76,77 through (1) CheckM, which uses Prodigal for gene calling, and (2) the more up to date HMMER3.2.2 on our gene calls described above. The latter generally produced slightly higher completeness and redundancy values (Supplementary Tables 8 and 9). For the expanded set of Asgard archaea genomes used for the phylogenomic analyses shown in Extended Data Fig. 4b, we applied the following filtering criteria: ≤100 contigs, >96% marker completeness and 20% sequence identity, >85% sequence alignment and 30% sequence identity, >90% sequence alignment and More

The prevalence of S. aureus isolation from the different examined samplesStaphylococcal infections represent a public health issue in hospitals and health care settings as well as a major economical and welfare problem in dairy animal farming25. The prevalence of S. aureus isolation from the farm under the study (Table 2) showed that 63 (33.1%) out of 190 different samples were bacteriologically positive. Moreover, the isolation was mainly obtained from manager swabs followed by milk machine swabs, nasal swabs and hand swabs (60.0, 53.3, 40.0 and 28.0%, respectively), and to a lesser extent in milk samples (24.0%). Meanwhile it was not isolated at any percent from water trough swabs, at X2 = 48.8, P  More

Generation timeDuring the infectious period, an infected individual may produce a secondary infection. However, the individual’s infectiousness is not constant during the infectious period, but it can be approximated by the probability distribution of the generation time (GT), which accounts for the time between the infection of a primary case and the infection of a secondary case. Unfortunately, such distribution is not as easy to estimate as that of the serial interval, which accounts for the time between the onset of symptoms in a primary case to the onset of symptoms of a secondary case. This is because the time of infection is more difficult to detect than the time of symptoms onset. Ganyani et al.27 developed a methodology to estimate the distribution of the GT from the distributions of the incubation period and the serial interval. Assuming an incubation period following a gamma distribution with a mean of 5.2 days and a standard deviation (SD) of 2.8 days, they estimated the serial interval from 91 and 135 pairs of documented infector-infectee in Singapore and Tianjin (China). Then, they found that the GT followed a gamma distribution with mean = 5.20 (95% CI = [3.78, 6.78]) days and SD = 1.72 (95% CI = [0.91, 3.93]) for Singapore (hereafter GT1), and with mean = 3.95 (95% CI = [3.01, 4.91]) days and SD = 1.51 (95% CI = [0.74, 2.97]) for Tianjin (hereafter GT2). Ng et al.28 applied the same methodology to 209 pairs of infector-infectee in Singapore and determined a gamma distribution with mean = 3.44 (95% CI = [2.79, 4.11]) days and SD 2.39 (95% CI = [1.27, 3.45]; hereafter GT3). Figure 3 shows the probability density functions (PDF) of such distributions, fGT. The differences between them are remarkable. For example, the 54.5%, 81.0%, and 80.7% of the contagions are produced in a pre-symptomatic stage (in the first 5.2 days after primary infection) assuming GT1, GT2, and GT3, respectively.Figure 3Probability density function of the generation time distribution, fGT(t), of GT1 (blue line; Singapore27), GT2 (yellow line; Tianjin27), GT3 (red line; Singapore28), and GTth (black line; theoretical distribution). Bars are the discretized version, (widetilde{{f_{GT} }}left( n right)), of the PDF of GTth.Full size imageTheoretically, assuming that the incubation periods of two individuals are independent and identically distributed, which is quite plausible, the expected/mean values of the GT and the serial interval should be equal29,30. The mean of the serial interval is easier to estimate than that of the GT. For that reason, we assume a mean serial interval as estimated from a meta-analysis of 13 studies involving a total of 964 pairs of infector-infectee, which is 4.99 days (95% CI = [4.17, 5.82])31, is more reliable than the aforementioned means of the GT. This value is within the error estimates of the means of GT1 and GT2, but not for GT3. Then, we construct a theoretical distribution for the GT that follows a gamma distribution (hereafter GTth) with mean = 4.99 days and SD = 1.88 days. This theoretical distribution can be seen in Fig. 3 and approximates the average PDF of three gamma distributions with mean = 4.99 and the SD of GT1, GT2, and GT3. We assume a conservative CI = [1.51, 2.39] for the theoretical SD, defined with the minimum and maximum SD values of GT1, GT2, and GT3. GTth shows 63.1% of pre-symptomatic contagions.
R

0

from r
In theory, the basic reproduction number R0 can be estimated as far as the intrinsic growth rate r, and the distributions of both the latent and infectious periods are known26,32,33,34. The latent period accounts for the period during which an infected individual cannot infect other individuals. It is observed in diseases for which the infectious period starts around the end of the incubation period, as happened with influenza35 and SARS36. However, from Fig. 3 it is inferred that COVID-19 is transmissible from the moment of infection, and we will assume a null latent period. Then, if the GT follows a gamma distribution, R0 can be estimated from the formulation of Anderson and Watson32, which was adapted to null latent periods by Yan26 as$$ R_{0} = frac{{mean_{GT} }}{{1 – left( {1 + mean_{GT} cdot r cdot frac{1}{{shape_{GT} }}} right)^{{ – shape_{GT} }} }} cdot r, $$
(4)
where meanGT is the mean GT and shapeGT is one of the two parameters defining the gamma distribution, which can be estimated as$$ shape_{GT} = frac{{left( {mean_{GT} } right)^{2} }}{{left( {SD_{GT} } right)^{2} }}. $$
(5)
For GTth, we get R0 = 1.50 (CI = [1.41, 1.61]) for REMEDID I(n) and R0 = 1.76 (CI = [1.60, 1.94]) for official I(n). For the other three GT distributions, R0 ranges from 1.39 (CI = [1.27, 1.58]) to 1.51 (CI = [1.34, 1.80]) for REMEDID I(n) and from 1.59 (CI = [1.40, 1.88]) to 1.78 (CI = [1.51, 2.23]) for official I(n) (Table 1). In all cases, R0 from GTth are within those from the three known GT distributions and indistinguishable from them within the error estimates. The lower (upper) bound of the CI is estimated as the minimum (maximum) R0 obtained from all the possible combinations of 100 evenly spaced values covering the CI of r, meanGT and SDGT. Then, following the Bonferroni correction, the reported CI present at least a 85% of confidence level for GT1, GT2, and GT3, but it cannot be assured for GTth since the CI of its SD is unknown. In general, all these R0 estimates are lower than those summarised by Park et al.20.Table 1 R0 and HIT values of the ancestral SARS-CoV-2 variant estimated from GT1, GT2, GT3, and GTth, and REMEDID and official infections. For date0, “Dec.” means December 2019, and “Jan.” means January 2020.Full size tableAlternatively, R0 can be estimated by applying the Euler–Lotka equation29,33,$$ R_{0} = frac{1}{{mathop smallint nolimits_{0}^{ + infty } e^{ – rt} cdot f_{GT} left( t right)dt}}. $$
(6)
In this case, we get values closer to previous estimates20. In particular, for GTth, we get R0 = 2.12 (CI = [1.81, 2.48]) for REMEDID I(n) and R0 = 2.92 (CI = [2.28, 3.75]) for official I(n). For the other three GT distributions, R0 ranges from 1.63 (CI = [1.43, 1.90]) to 2.21 (CI = [1.59, 2.95]) for REMEDID I(n) and from 1.97 (CI = [1.59, 2.54]) to 3.11 (CI = [1.84, 4.90]) for official I(n) (Table 1). The CI are estimated as in Eq. (4).R0 from a dynamical modelWe designed a dynamic model with Susceptible-Infected-Recovered (SIR) as stocks that accounts for the infectiousness of the infectors. Such a model is a generalisation of the Susceptible-Exposed-Infected-Recovered (SEIR) model37. Births, deaths, immigration and emigration are ignored, which seems reasonable since the timescale of the outbreak is too short to produce significant demographic changes. For the sake of simplicity, the recovered stock includes recoveries and fatalities, and it is denoted as R(t). A random mixing population is assumed, that is a population where contacts between any two people are equally probable. Time is discretized in days, so the real time variable t is replaced by the integer variable n. As a consequence, the derivatives in the differential equations defining the dynamic model explained below are discrete derivatives.The size of the population is fixed at N = 100,000, and then, for any day n we get$$ tilde{S}left( n right) + left( {mathop sum limits_{k = 0}^{20} tilde{I}left( n-k right)} right) + tilde{R}left( n right) = N, $$
(7)
where (tilde{S}left( n right)), (tilde{I}left( n right)), and (tilde{R}left( n right)) are the discretized versions of S(t), I(t), and R(t) and (tilde{I}) is assumed to be null for negative integers. The summation is a consequence of the infectiousness, which is approximated according to the GT, whose PDF is discretized as$$ widetilde{{f_{GT} }}left( n right) = mathop smallint limits_{n – 1}^{n} f_{GT} left( t right) dt, $$
(8)
from n = 1 to 20. Figure 3 shows (widetilde{{f_{GT} }}left( n right)) for GTth. Truncating at n = 20 accounts for 99.99% of the area below the PDF of all the GT. Then, an infected individual at day n0 is expected to produce on average$$ widetilde{{R_{e} }}left( {n_{0} + n} right) cdot widetilde{{f_{GT} }}left( n right) $$
(9)
infections n days later, where (widetilde{{R_{e} }}left( n right)) is the discretized version of Re(t). From this expression, it is obvious that values of (widetilde{{R_{e} }}left( n right) < 1) will produce a decline of infections. Conversely, infections at day n0 are produced by all individuals infected during the previous 20 days as$$ tilde{I}(n_{0} ) = tilde{R}_{e} left( {n_{0} } right) cdot left( {mathop sum limits_{n = 1}^{20} tilde{I}left( {n_{0} - n} right) cdot widetilde{{f_{GT} }}left( n right)} right), $$ (10) whose continuous version has been reported in previous studies29,38. The expression in brackets is called total infectiousness of infected individuals at day n039. According to Eq. (1), Eq. (10) can be expressed in terms of R0 as$$ tilde{I}(n_{0} ) = R_{0} cdot frac{{tilde{S}left( {n_{0} } right)}}{N} cdot left( {mathop sum limits_{n = 1}^{20} tilde{I}left( {n_{0} - n} right) cdot widetilde{{f_{GT} }}left( n right)} right). $$ (11) As we want a dynamic model capable of providing (tilde{I}left( {n_{0} } right)) from the stocks at time step n0 − 1, we replaced (tilde{S}left( {n_{0} } right)) by (tilde{S}left( {n_{0} - 1} right)) in Eq. (11). This assumption makes sense in a discrete domain since the infections at time n0 take place in the susceptible population at time n0 − 1. Then, assuming that all stocks are set to zero for negative integers, our dynamic model can be expressed in terms of Eq. (7) and the following differential equations:$$ delta tilde{I}(n_{0} ) = R_{0} cdot frac{{tilde{S}left( {n_{0} - 1} right)}}{N} cdot left( {mathop sum limits_{n = 1}^{20} tilde{I}left( {n_{0} - n} right) cdot widetilde{{{text{f}}_{GT} }}left( n right)} right) - tilde{I}(n_{0} - 1), $$ (12) $$ delta tilde{S}left( {n_{0} } right) = {-}tilde{I}left( {n_{0} } right), $$ (13) $$ delta tilde{R}left( {n_{0} } right) = tilde{I}left( {n_{0} - 21} right), $$ (14) where (delta tilde{I}), (delta tilde{S}), and (delta tilde{R}) are the (discrete) derivatives of (tilde{I}), (tilde{S}), and (tilde{R}), respectively. Applying the initial conditions (tilde{S}left( 0 right) = N - 1), (tilde{I}left( 0 right) = 1), and (tilde{R}left( 0 right) = 0), it is assumed that the outbreak was produced by only one infector. The latter is not true in Spain, since several independent introductions of SARS-CoV-2 were detected40. However, for modelling purposes it is equivalent to introducing a single infection at day 0 or M infections produced by the single infection n days later. Then, the date of the initial time n = 0 is accounted as a parameter date0, which is optimised, as well as R0, to minimise the root-mean square of the residual between the model simulated (tilde{I}left( n right)) and the REMEDID and official I(n) for the period from date0 to n0.The model was implemented in Stella Architect software v2.1.1 (www.iseesystems.com) and exported to R software v4.1.1 with the help of deSolve (v1.28) and stats (v4.1.1) packages, and the Brent optimisation algorithm was implemented. For REMEDID I(n) and GTth, we obtained date0 = 13 December 2019 and R0 = 2.71 (CI = [2.33, 3.15]). Optimal solutions combine lower/higher R0 and earlier/later date0 (Fig. 4), which highlights the importance of providing an accurate first infection date to estimate R0. When the other three GT distributions were considered, we obtained similar date0, ranging from 12 to 17 December 2019, and R0 values ranging from 2.08 (CI = [1.86, 2.42]) to 2.85 (CI = [2.05, 3.25]; see Table 1). For official infections, date0 was set to 1 January 2020 for all cases, and R0 ranged from 1.81 (CI = [1.64, 2.07]) to 2.41 (CI = [1.80, 2.91]). The CI are estimated as in Eq. (4).Figure 4Root-mean square (RMS) of the residuals between infections from the model, which depends on date0 (x-axis) and R0 (y-axis), and REMEDID (from MoMo ED) and official infections. Parameters optimizing the model are highlighted in purple. RMS larger than 1275 (left panel) and 103 (right panel) are saturated in white.Full size image More

1.Federal Office for the Environment (FOEN). Action Plan for the Swiss Biodiversity Strategy (FOEN, Bern, 2017).2.Geldmann, J. & González-Varo, J. P. Conserving honey bees does not help wildlife. Science 359, 392–393 (2018).Article 

Google Scholar 
3.Egerer, M. & Kowarik, I. Confronting the modern gordian knot of urban beekeeping. Trends Ecol. Evol. 35, 956–959 (2020).Article 

Google Scholar 
4.Herrera, C. M. Gradual replacement of wild bees by honeybees in flowers of the Mediterranean Basin over the last 50 years. Proc. R. Soc. B Biol. Sci. 287, 16–20 (2020).
Google Scholar 
5.Torné-Noguera, A., Rodrigo, A., Osorio, S. & Bosch, J. Collateral effects of beekeeping: impacts on pollen-nectar resources and wild bee communities. Basic Appl. Ecol. 17, 199–209 (2016).Article 

Google Scholar 
6.Magrach, A., González-Varo, J. P., Boiffier, M., Vilà, M. & Bartomeus, I. Honeybee spillover reshuffles pollinator diets and affects plant reproductive success. Nat. Ecol. Evol. 1, 1299–1307 (2017).Article 

Google Scholar 
7.Prendergast, K. S., DIxon, K. W. & Bateman, P. W. Interactions between the introduced European honey bee and native bees in urban areas varies by year, habitat type and native bee guild. Biol. J. Linn. Soc. 133, 725–743 (2021).Article 

Google Scholar 
8.Herbertsson, L., Lindström, S. A. M., Rundlöf, M., Bommarco, R. & Smith, H. G. Competition between managed honeybees and wild bumblebees depends on landscape context. Basic Appl. Ecol. 17, 609–616 (2016).Article 

Google Scholar 
9.Stevenson, P. C. et al. The state of the world’s urban ecosystems: what can we learn from trees, fungi, and bees? Plants People Planet 2, 482–498 (2020).Article 

Google Scholar 
10.Ropars, L., Dajoz, I., Fontaine, C., Muratet, A. & Geslin, B. Wild pollinator activity negatively related to honey bee colony densities in urban context. PLoS ONE 14, e0222316 (2019).CAS 
Article 

Google Scholar 
11.Venter, Z. S. & Sydenham, M. A. K. Continental-scale land cover mapping at 10 m resolution over Europe (ELC10). Remote Sens. 13, 2301 (2021).Article 

Google Scholar 
12.Hardin, G. The tragedy of the commons. Science 162, 1243–1248 (1968).CAS 
Article 

Google Scholar 
13.Tew, N. E. et al. Quantifying nectar production by flowering plants in urban and rural landscapes. J. Ecol. 109, 1747–1757 (2021).Article 

Google Scholar 
14.Baldock, K. C. R. et al. A systems approach reveals urban pollinator hotspots and conservation opportunities. Nat. Ecol. Evol. 3, 363–373 (2019).Article 

Google Scholar 
15.Casanelles-Abella, J. et al. Applying predictive models to study the ecological properties of urban ecosystems: A case study in Zürich, Switzerland. Landsc. Urban Plan. 214, 104137 (2021).Article 

Google Scholar 
16.IPBES. Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES, 2019).17.IUCN. IUCN’s Key Messages. First Draft of the Post-2020 Global Biodiversity Framework. In Convention on Biological Diversity Third meeting of the Open-Ended Working Group on the Post-2020 Global Biodiversity Framework (OEWG3) 15 (IUCN, 2021).18.Baldock, K. C. R. Opportunities and threats for pollinator conservation in global towns and cities. Curr. Opin. Insect Sci. 38, 63–71 (2020).Article 

Google Scholar 
19.Henry, M. & Rodet, G. The apiary influence range: a new paradigm for managing the cohabitation of honey bees and wild bee communities. Acta Oecologica 105, 103555 (2020).Article 

Google Scholar 
20.Ignatieva, M. & Hedblom, M. An alternative urban green carpet. Science 362, 148–149 (2018).CAS 
Article 

Google Scholar 
21.Vega, K. A. & Küffer, C. Promoting wildflower biodiversity in dense and green cities: The important role of small vegetation patches. Urban For. Urban Green. 62, 127165 (2021).Article 

Google Scholar 
22.Fabián, D., González, E., Sánchez Domínguez, M. V., Salvo, A. & Fenoglio, M. S. Towards the design of biodiverse green roofs in Argentina: assessing key elements for different functional groups of arthropods. Urban For. Urban Green. 61 (2021).23.Vega, K. A., Schläpfer-Miller, J. & Kueffer, C. Discovering the wild side of urban plants through public engagement. Plants People Planet 3, 389–401 (2021).Article 

Google Scholar 
24.Von Büren, R. S., Oehen, B., Kuhn, N. J. & Erler, S. High-resolution maps of Swiss apiaries and their applicability to study spatial distribution of bacterial honey bee brood diseases. PeerJ 2019, 1–21 (2019).
Google Scholar 
25.QGIS Development Team. QGIS Geographic Information System (Open Source Geospatial Foundation Project, 2020).26.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, 2019).27.R Studio Team. R studio: Integrated Development for R (RStudio, Boston, 2020).28.Casanelles-Abella, J. & Moretti M. Challenging the sustainability of urban beekeeping: evidence from Swiss cities. Envidat. https://doi.org/10.16904/envidat.239 (2021).29.Casanelles-Abella, J. Code for the paper: Challenging the sustainability of urban beekeeping: evidence from Swiss cities (v. 1.0). Zenodo https://doi.org/10.5281/zenodo.5618254 (2021).30.Swiss Fedearl Office of Topography Swisstopo. SWISSIMAGE 25. https://www.swisstopo.admin.ch/en/geodata/images/ortho/swissimage25.html#links (Swisstopo, 2021). More

1.Lindow SE, Brandl MT. Microbiology of the phyllosphere. Appl Environ Microbiol. 2003;69:1875LP–1883.
Google Scholar 
2.Morella NM, Zhang X, Koskella B. Tomato seed-associated bacteria confer protection of seedlings against foliar disease caused by Pseudomonas syringae. Phytobiomes J. 2019;3:177–90.
Google Scholar 
3.Innerebner G, Knief C, Vorholt JA. Protection of Arabidopsis thaliana against leaf-pathogenic Pseudomonas syringae by Sphingomonas strains in a controlled model system. Appl Environ Microbiol. 2011;77:3202–10.CAS 
PubMed 
PubMed Central 

Google Scholar 
4.Fu S-F, Sun P-F, Lu H-Y, Wei J-Y, Xiao H-S, Fang W-T, et al. Plant growth-promoting traits of yeasts isolated from the phyllosphere and rhizosphere of Drosera spatulata Lab. Fungal Biol. 2016;120:433–48.CAS 
PubMed 

Google Scholar 
5.Laforest-Lapointe I, Paquette A, Messier C, Kembel SW. Leaf bacterial diversity mediates plant diversity and ecosystem function relationships. Nature. 2017;546:145–7.CAS 
PubMed 

Google Scholar 
6.Lindow SE, Leveau JHJ. Phyllosphere microbiology. Curr Opin Biotechnol. 2002;13:238–43.CAS 
PubMed 

Google Scholar 
7.Fürnkranz M, Wanek W, Richter A, Abell G, Rasche F, Sessitsch A. Nitrogen fixation by phyllosphere bacteria associated with higher plants and their colonizing epiphytes of a tropical lowland rainforest of Costa Rica. ISME J. 2008;2:561–70.PubMed 

Google Scholar 
8.Ottesen AR, Gorham S, Reed E, Newell MJ, Ramachandran P, Canida T, et al. Using a control to better understand phyllosphere microbiota. PLoS ONE. 2016;11:e0163482.PubMed 
PubMed Central 

Google Scholar 
9.Jones JDG, Dangl JL. The plant immune system. Nature. 2006;444:323–9.CAS 
PubMed 

Google Scholar 
10.Bodenhausen N, Bortfeld-miller M, Ackermann M, Vorholt JA. A synthetic community approach reveals plant genotypes affecting the phyllosphere microbiota. PLoS Biol. 2014; 10. https://doi.org/10.1371/journal.pgen.1004283.11.Horton MW, Bodenhausen N, Beilsmith K, Meng D, Muegge BD, Subramanian S, et al. Genome-wide association study of Arabidopsis thaliana leaf microbial community. Nat Commun. 2014;5:5320.PubMed 

Google Scholar 
12.Hacquard S, Spaepen S, Garrido-Oter R, Schulze-Lefert P. Interplay between innate immunity and the plant microbiota. Annu Rev Phytopathol. 2017;55:565–89.CAS 
PubMed 

Google Scholar 
13.Zhalnina K, Louie KB, Hao Z, Mansoori N, Nunes U, Shi S et al. Dynamic root exudate chemistry and substrate preferences drive patterns in rhizosphere microbial community assembly. Nat Microbiol. 2018. https://doi.org/10.1038/s41564-018-0129-3.14.Humphrey PT, Whiteman NK. Insect herbivory reshapes a native leaf microbiome. Nat Ecol Evol. 2020;4:221–9.PubMed 
PubMed Central 

Google Scholar 
15.Yadav RKP, Karamanoli K, Vokou D. Bacterial colonization of the phyllosphere of Mediterranean perennial species as influenced by leaf structural and chemical features. Micro Ecol. 2005;50:185–96.CAS 

Google Scholar 
16.Morella NM, Weng FCH, Joubert PM, Metcalf CJE, Lindow S, Koskella B. Successive passaging of a plant-associated microbiome reveals robust habitat and host genotype-dependent selection. Proc Natl Acad Sci USA. 2020;117:1148–59.CAS 
PubMed 

Google Scholar 
17.Wagner MR, Busby PE, Balint-Kurti P. Analysis of leaf microbiome composition of near-isogenic maize lines differing in broad-spectrum disease resistance. N Phytol. 2019;225:2152–65.
Google Scholar 
18.Wagner MR, Lundberg DS, Del Rio TG, Tringe SG, Dangl JL, Mitchell-Olds T. Host genotype and age shape the leaf and root microbiomes of a wild perennial plant. Nat Commun. 2016;7:1–15.CAS 

Google Scholar 
19.Horner-Devine MC, Bohannan BJM. Phylogenetic clustering and overdispersion in bacterial communities. Ecology. 2006;87:S100–8.PubMed 

Google Scholar 
20.Kembel SW, O’Connor TK, Arnold HK, Hubbell SP, Wright SJ, Green JL. Relationships between phyllosphere bacterial communities and plant functional traits in a neotropical forest. Proc Natl Acad Sci USA 2014; 1–6.21.Burns AR, Stephens WZ, Stagaman K, Wong S, Rawls JF, Guillemin K, et al. Contribution of neutral processes to the assembly of gut microbial communities in the zebrafish over host development. ISME J. 2016;10:655–64.CAS 
PubMed 

Google Scholar 
22.Sloan WT, Woodcock S, Lunn M, Head IM, Curtis TP. Modeling taxa-abundance distributions in microbial communities using environmental sequence data. Micro Ecol. 2007;53:443–55.
Google Scholar 
23.Laforest-Lapointe I, Messier C, Kembel SW. Host species identity, site and time drive temperate tree phyllosphere bacterial community structure. Microbiome 2016; 1–10.24.Schlaeppi K, Dombrowski N, Oter RG, Ver Loren van Themaat E, Schulze-Lefert P. Quantitative divergence of the bacterial root microbiota in Arabidopsis thaliana relatives. Proc Natl Acad Sci USA. 2014;111:585LP–592.
Google Scholar 
25.Gallart M, Adair KL, Love J, Meason DF, Clinton PW, Xue J, et al. Host genotype and nitrogen form shape the root microbiome of Pinus radiata. Micro Ecol. 2018;75:419–33.CAS 

Google Scholar 
26.Hambäck PA, Inouye BD, Andersson P, Underwood N. Effects of plant neighborhoods on plant–herbivore interactions: resource dilution and associational effects. Ecology. 2014;95:1370–83.PubMed 

Google Scholar 
27.Underwood N, Inouye BD, Hambäck PA. A conceptual framework for associational effects: when do neighbors matter and how would we know? Q Rev Biol. 2014;89:1–19.PubMed 

Google Scholar 
28.Barbosa P, Hines J, Kaplan I, Martinson H, Szczepaniec A, Szendrei Z. Associational resistance and associational susceptibility: having right or wrong neighbors. Annu Rev Ecol Evol Syst. 2009;40:1–20.
Google Scholar 
29.Janzen DH. Herbivores and the number of tree species in tropical forests. Am Nat. 1970;104:501–28.
Google Scholar 
30.Connell JH. On the role of natural enemies in preventing competitive exclusion in some marine animals and in rain forest trees. in Den Boer PJ, Gradwell G, editors. Dynamics of populations. PUDOC, 1971, p. 298–312.31.Mangan SA, Schnitzer SA, Herre EA, Mack KML, Valencia MC, Sanchez EI, et al. Negative plant–soil feedback predicts tree-species relative abundance in a tropical forest. Nature. 2010;466:752–5.CAS 
PubMed 

Google Scholar 
32.Miller EC, Perron GG, Collins CD. Plant‐driven changes in soil microbial communities influence seed germination through negative feedbacks. Ecol Evol. 2019;0:1–14.
Google Scholar 
33.Antonovics J, Ellstrand NC. Experimental studies of the evolutionary significance of sexual reproduction. I. A test of the frequency-dependent selection hypothesis. Evolution. 1984;38:103–15.PubMed 

Google Scholar 
34.Ellstrand NC, Antonovics J. Experimental studies of the evolutionary significance of sexual reproduction II. A test of the density-dependent selection hypothesis. Evolution. 1985;39:657–66.PubMed 

Google Scholar 
35.Naeem S, Tjossem SF, Byers D, Bristow C, Li S. Plant neighborhood diversity and production. Ecoscience. 1999;6:355–65.
Google Scholar 
36.Worrich A, Musat N. Associational effects in the microbial neighborhood. ISME J 2019; 2143–9.37.Copeland JK, Yuan L, Layeghifard M, Wang PW, Guttman DS. Seasonal community succession of the phyllosphere microbiome. Mol Plant Microbe Interact. 2015;28:274–85.CAS 
PubMed 

Google Scholar 
38.Lajoie G, Kembel SW. Host neighborhood shapes bacterial community assembly and specialization on tree species across a latitudinal gradient. Ecol Monogr. 2021;0:1–18.
Google Scholar 
39.Lymperopoulou D, Adams R, Lindow SE, Löffler F. Contribution of vegetation to the microbial composition of nearby outdoor air. Appl Environ Microbiol. 2016;82:3822–33.CAS 
PubMed 
PubMed Central 

Google Scholar 
40.Lindow SE, Andersen G. Influence of immigration on epiphytic bacterial populations on navel orange leaves. Appl Environ Microbiol. 1996;62:2978–87.CAS 
PubMed 
PubMed Central 

Google Scholar 
41.Massoni J, Bortfeld-miller M, Widmer A, Vorholt JA. Capacity of soil bacteria to reach the phyllosphere and convergence of floral communities despite soil microbiota variation. Proc Natl Acad Sci USA 2021; 118. https://doi.org/10.1073/pnas.2100150118.42.Leibold MA, Holyoak M, Mouquet N, Amarasekare P, Chase JM, Hoopes MF, et al. The metacommunity concept: a framework for multi-scale community ecology. Ecol Lett. 2004;7:601–13.
Google Scholar 
43.Fodelianakis S, Lorz A, Valenzuela-cuevas A, Barozzi A, Booth JM, Daffonchio D. Dispersal homogenizes communities via immigration even at low rates in a simplified synthetic bacterial metacommunity. Nat Commun. 2019;10:1–12.
Google Scholar 
44.Burns AR, Miller E, Agarwal M, Rolig AS, Milligan-Myhre K, Seredick S et al. Interhost dispersal alters microbiome assembly and can overwhelm host innate immunity in an experimental zebrafish model. Proc Natl Acad Sci USA 2017;114. https://doi.org/10.1073/pnas.1702511114.45.Chelius MK, Triplett EW. The diversity of archaea and bacteria in association with the roots of Zea mays L. Micro Ecol. 2001;41:252–63.CAS 

Google Scholar 
46.Bodenhausen N, Horton MW, Bergelson J. Bacterial communities associated with the leaves and the roots of Arabidopsis thaliana. PLoS ONE. 2013;8:e56329.CAS 
PubMed 
PubMed Central 

Google Scholar 
47.Lundberg DS, Yourstone S, Mieczkowski P, Jones CD, Dangl JL. Practical innovations for high-throughput amplicon sequencing. Nat Methods. 2013;10:999–1002.CAS 
PubMed 

Google Scholar 
48.Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581.CAS 
PubMed 
PubMed Central 

Google Scholar 
49.R Core Team. R: A language and environment for statistical computing. 2020. http://cran.r-project.org.50.Morgan M, Anders S, Lawrence M, Aboyoun P, Pagès H, Gentleman R. ShortRead: a bioconductor package for input, quality assessment and exploration of high-throughput sequence data. Bioinformatics. 2009;25:2607–8.CAS 
PubMed 
PubMed Central 

Google Scholar 
51.Pagès H, Aboyoun P, Gentleman R, DebRoy S. Biostrings: efficient manipulation of biological strings. 2020.52.McMurdie P, Holmes S. phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 2013; 8. https://doi.org/10.1371/journal.pone.0061217tle.53.Wang Q, Garrity GM, Tiedje JM, Cole JR. Naive bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microbiol. 2007;73:5261–7.CAS 
PubMed 
PubMed Central 

Google Scholar 
54.Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:D590–D596.CAS 
PubMed 

Google Scholar 
55.Davis NM, Proctor DM, Holmes SP, Relman DA, Callahan BJ. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome. 2018;6:226.PubMed 
PubMed Central 

Google Scholar 
56.Morella NM, Gomez AL, Wang G, Leung MS, Koskella B. The impact of bacteriophages on phyllosphere bacterial abundance and composition. Mol Ecol. 2018;27:2025–38.PubMed 

Google Scholar 
57.Oksanen J, Blanchet FG, Roeland K, Legendre P, Minchin P, O’Hara RB et al. vegan: Community ecology package. 2015. http://cran.r-project.org.58.Anderson MJ. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 2001;26:32–46.
Google Scholar 
59.De Cáceres M, Legendre P. Associations between species and groups of sites: indices and statistical inference. Ecology. 2009;90:3566–74.PubMed 

Google Scholar 
60.Ersts PJ. Geographic Distance Matrix Generator. http://biodiversityinformatics.amnh.org/open_source/gdmg.61.Sprockett D. reltools: Microbiome Amplicon Analysis and Visualization. 2021.62.Sloan WT, Lunn M, Woodcock S, Head IM, Nee S, Curtis TP. Quantifying the roles of immigration and chance in shaping prokaryote community structure. Environ Microbiol. 2006;8:732–40.PubMed 

Google Scholar 
63.Wright ES. Using DECIPHER v2.0 to analyze big biological sequence data in R. R J. 2016;8:352–9.
Google Scholar 
64.Schliep KP. phangorn: phylogenetic analysis in R. Bioinformatics. 2011;27:592–3.CAS 
PubMed 

Google Scholar 
65.Kembel SW, Cowan PD, Helmus MR, Cornwell WK, Morlon H, Ackerly DD, et al. Picante: R tools for integrating phylogenies and ecology. Bioinformatics. 2010;26:1463–4.CAS 
PubMed 

Google Scholar 
66.Koskella B. The phyllosphere. Curr Biol. 2020;30:R1143–R1146.CAS 
PubMed 

Google Scholar 
67.Chaparro JM, Badri DV, Vivanco JM. Rhizosphere microbiome assemblage is affected by plant development. ISME J. 2014;8:790–803.CAS 
PubMed 

Google Scholar 
68.İnceoğlu Ö, Al-Soud WA, Salles JF, Semenov AV, van Elsas JD. Comparative analysis of bacterial communities in a potato field as determined by pyrosequencing. PLoS ONE. 2011;6:e23321.PubMed 
PubMed Central 

Google Scholar 
69.Christian N, Herre EA, Mejia LC, Clay K. Exposure to the leaf litter microbiome of healthy adults protects seedlings from pathogen damage. Proc R Soc B Biol Sci. 2017;284:20170641.
Google Scholar 
70.Leigh EG, Davidar P, Dick CW, Terborgh J, Puyravaud J-P, ter Steege H, et al. Why do some tropical forests have so many species of trees? Biotropica. 2004;36:447–73.
Google Scholar 
71.Hyatt LA, Rosenberg MS, Howard TG, Bole G, Fang W, Anastasia J, et al. The distance dependence prediction of the Janzen-Connell hypothesis: a meta-analysis. Oikos. 2003;103:590–602.
Google Scholar 
72.Carson W, Anderson J, Leigh E, Schnitzer S. Challenges associated with testing and falsifying the Janzen_Connell hypothesis: a review and critique. In: Carson W, Schnitzer SA, editors. Tropical forest community ecology. Wiley Blackwell; 2008. p. 210–41. More

Genome sequencing and assemblyGenome assemblies ranged in size from 799.9 Mbp in Melanocharis versteri to 1053.5 Mbp in Sericornis nouhuysi. The number of scaffolds ranged from 14,086 scaffolds in Melipotes ater to 87,957 scaffolds in Ficedula hyperythra and N50 ranged between ca. 40 Kbp to and 25 Mbp. Benchmarking Universal Single-Copy Orthologs (BUSCO) analyses of genome completeness ranged from a high proportion of complete BUSCOs in Melipotes ater, 86.8% to only 66.7% complete BUSCOs in Rhipidura albolimbata. For most species, the proportions of complete BUSCOs were 75–80%. Overall, the proportion of missing BUSCOs was low, ranging from 6.6% in Melipotes ater to 15.2% in Rhipidura albolimbata (see Supplementary Table 1 for all genome assembly statistics and Supplementary Fig. 1 for the number of SNP variants per species).Kinship analyses of individuals within populationsSampling of closely related individuals can dramatically bias estimates of population structure and demographics. Two Pachycephala schlegelii individuals (A117 and A118) showed a pairwise kinship coefficient of 0.144, indicative of being half-siblings. The two individuals were collected at the same locality on the same date. Similarly, two Ifrita kowaldi individuals (D116 and D117) showed a pairwise kinship coefficient of 0.135, also suggestive of being half-siblings. In this case, the individuals were collected on the same sampling locality on two consecutive days. To not bias downstream demographic analyses, one of the P. schlegelii (A118) and one of the I. kowaldi (D117) individuals were excluded from all subsequent analyses. For all other species, no closely related individuals were identified.Genetic differentiationEstimated levels of differentiation between populations were initially based on three approaches; (i) calculation of FST (the fixation index), which quantifies the degree of genetic differentiation between populations based on the variation in allele frequencies, ranging between 0 (complete mixing of individuals) and 1 (complete separation of populations) (Fig. 1), (ii) Standardized pairwise FST used to conduct a Principal Component Analysis (PCA) in order to visualize population structure (Supplementary Fig. 1) and (iii) Admixture analysis as implemented in STRUCTURE (a clustering algorithm that infers the most likely number of groups [K]), in which individuals are grouped into clusters according to the proportion of their ancestry components (Supplementary Fig. 1). As a preliminary analysis, we calculated FST and constructed PCA plots for the four congeneric (incl. Sericornis/Aethomyias [until recently placed in the genus Sericornis]) species pairs in our study (Supplementary Fig. 2), which were aligned using the same reference genome. This was done to ascertain that no samples had been misidentified and to gauge levels of differentiation between distinct species. All species were genetically well separated and FST values ranged from 0.08 for the two Ptiloprora species to 0.20 for the two Ficedula species.For five out of six species from Buru/Seram, genetic differentiation (FST) was high between islands (Fig. 1), and comparable to differentiation between named congeneric species in this study (e.g. Ptiloprora and Melipotes); Ceyx lepidus (FST = 0.16), Thapsinillas affinis (FST = 0.15), Ficedula buruensis (FST = 0.13) and Pachycephala macrorhyncha (FST = 0.09). In contrast, differentiation in Ficedula hyperythra was consistent with population-level differentiation (FST = 0.04). In all cases, individuals from Buru and Seram were clearly differentiated in the PCA and STRUCTURE plots (Supplementary Fig. 1A). For Ceyx lepidus, Ficedula buruensis and Pachycephala macrorhyncha, samples were collected at multiple elevations and we therefore calculated genetic differentiation between elevations (Buru: 1097 m versus 1435 m and Seram: 1000 m versus 1300 m) to determine any potential parapatric differentiation along the gradients. In all possible comparisons, FST values did not differ significantly from 0. Moreover, PCA plots showed that samples did not cluster according to elevation (Supplementary Fig. 3A).Three of the thirteen New Guinean population pairs occurring in Mount Wilhelm and Huon showed relatively high genetic divergences: Melipotes fumigatus/ater (FST = 0.08), Paramythia montium (FST = 0.09) and Ifrita kowaldi (FST = 0.07) (Fig. 1) with populations clearly separated (Supplementary Fig. 1). By contrast, the two lowland species Toxorhamphus novaeguineae and Melilestes megarhynchus showed little genetic differentiation, FST = 0.00. For the remaining species, genetic differentiation between Mount Wilhelm and Huon ranged between FST = 0.01–0.05. Despite this moderate level of genetic differentiation, the populations of Mount Wilhelm and Huon could be clearly distinguished in the PCA plots. In all cases STRUCTURE suggested a scenario with K = 2 with some mixing of individuals, except for Rhipidura albolimbata, in which K = 1 was suggested.For five bird species we included an additional population from Mount Scratchley, which is also situated in the Central Range but ~400 km to the southeast of Mount Wilhelm. Genetic differentiation of this population from the other two populations was comparable with that between Mount Wilhelm and Huon. The highest genetic differentiation was found in Paramythia montium (FST = 0.10 both between Mount Wilhelm and Mount Scratchley and between Huon and Mount Scratchley). In the case of Peneothello sigillata, the Mount Scratchley population appeared genetically well-differentiated from both the populations of Mount Wilhelm (FST = 0.06) and Huon (FST = 0.07). In both cases, STRUCTURE suggested a scenario of K = 3, with individual assignments matching the three geographically circumscribed populations. For Pachycephala schlegelii, genetic differentiation was relatively high between Huon and Mount Scratchley (FST = 0.05), but low between Mount Wilhelm and Mount Scratchley (FST = 0.01). Accordingly, STRUCTURE suggested a scenario with K = 2 groups. For the remaining two species Sericornis nouhuysi showed some differentiation (FST = 0.03) between Mount Wilhelm and Huon and Aethomyias papuensis showed minor differentiation (FST = 0.02 between Mount Scratchley and Huon (Supplementary Table 2), but for both species, STRUCTURE suggested a scenario of K = 2 with considerable mixing of individuals between populations.Samples from Mount Wilhelm were collected at elevations ranging from 1700 to 3700 m, again allowing us to test for differences within populations on a single slope, a finding that would be consistent with incipient parapatric speciation. No species showed significant differences in FST when comparing individuals from different elevations, and concordantly there was little clustering of individuals by elevation in the PCA plots. Even when individuals were collected as far as 2000 elevational meters apart (as in the case of Origma robusta), genetic differentiation was low (FST = 0.01). In Huon, all samples were collected at the same elevation, except for Ifrita kowaldi, for which genetic differentiation of FST = 0.03 was found between individuals collected at 2300 m and 2950 m (Supplementary Fig. 3B, Supplementary Table 2). These analyses however, suffer from very small sample sizes that hinder a thorough analysis of parapatric speciation events. Furthermore, we note that divergence with gene flow may not manifest as a genome-wide phenomenon (at least, not until the taxa are so differentiated that gene flow has ceased). Instead, it may proceed via selection acting to create small ‘islands of differentiation’ within the genome against a background of negligible differentiation22,23. Such analyses require large sample sizes and are therefore not possible herein.Correlations between genetic divergence and elevationIf lineages colonize mountains from the lowlands, followed by range contraction and differentiation in the highlands, we would expect a signature of larger genetic differentiation (FST) between populations inhabiting higher elevations. We found no relationship between genetic differentiation (FST) and the altitudinal floor (the lowest elevation at which a species/population occurs) for the five Moluccan species, but for all New Guinean taxa with the exception of Melipotes fumigatus/ater we found a significant positive correlation (r = 0.83, p  More

1.Kawecki, T. J. & Ebert, D. Conceptual issues in local adaptation. Ecol. Lett. 7, 1225–1241 (2004).
Google Scholar 
2.Fisher, T. W. et al. Handbook of Biological Control: Principles and Applications of Biological Control (Academic Press, London, 1999).
Google Scholar 
3.Dyck, V. A., Hendrichs, J. & Robinson, A. S. Sterile insect technique: Principles and practice in area-wide integrated pest management. In Sterile Insect Technique: Principles and Practice in Area-Wide Integrated Pest Management. https://doi.org/10.1007/1-4020-4051-2. (2005)4.Etges, W. J. & Noor, M. A. F. Genetics of Mate Choice: From Sexual Selection to Sexual Isolation. (Kluwer Academic Publishers, 2002).5.Harbach, R. E. Review of the internal classification of the genus Anopheles (Diptera: Culicidae): The foundation for comparative systematics and phylogenetic research. Bull. Entomol. Res. 84, 331–342 (1994).
Google Scholar 
6.Rogers, D. J., Randolph, S. E., Snow, R. W. & Hay, S. I. Satellite imagery in the study and forecast of malaria. Nature 415, 710–715 (2002).CAS 
PubMed 
PubMed Central 

Google Scholar 
7.Bayoh, M. N. N., Thomas, C. J. J. & Lindsay, S. W. W. Mapping distributions of chromosomal forms of Anopheles gambiae in West Africa using climate data. Med. Vet. Entomol. 15, 267–274 (2001).CAS 
PubMed 

Google Scholar 
8.Namountougou, M. et al. Multiple insecticide resistance in Anopheles gambiae s. l. Populations from Burkina Faso. West Africa. PLoS One 7, e48412 (2012).CAS 
PubMed 
ADS 

Google Scholar 
9.Benedict, M. Q. & Robinson, A. S. The first releases of transgenic mosquitoes: An argument for the sterile insect technique. Trends Parasitol. 19, 349–355 (2003).PubMed 

Google Scholar 
10.Maïga, H. et al. Mating competitiveness of sterile male Anopheles coluzzii in large cages. Malar. J. 13, 460 (2014).PubMed 
PubMed Central 

Google Scholar 
11.Clements, A. N. The Biology of Mosquitoes. Sensory Reception and Behaviour Behaviour, Vol. 2. (Wallingford, 1999).12.Doug, P. et al. Genetic and environmental factors associated with laboratory rearing affect survival and assortative mating but not overall mating success in Anopheles gambiae Sensu Stricto. PLoS One 8, e82631 (2013).
Google Scholar 
13.Baeshen, R. et al. Differential effects of inbreeding and selection on male reproductive phenotype associated with the colonization and laboratory maintenance of Anopheles gambiae. Malar. J. 13, 19 (2014).PubMed 
PubMed Central 

Google Scholar 
14.Bargielowski, I., Kaufmann, C., Alphey, L., Reiter, P. & Koella, J. Flight performance and teneral energy reserves of two genetically-modified and one wild-type strain of the yellow fever mosquito Aedes aegypti. Vector-Borne Zoonotic Dis. 12, 1053–1058 (2012).PubMed 
PubMed Central 

Google Scholar 
15.Harris, A. F. et al. Field performance of engineered male mosquitoes. Nat. Biotechnol. 29, 1034–1037 (2011).CAS 
PubMed 

Google Scholar 
16.Alphey, L. et al. Genetic control of Aedes mosquitoes. Pathogens and Global Health 107, 170–179 (2013).PubMed 
PubMed Central 

Google Scholar 
17.Lehmann, T. et al. Tracing the origin of the early wet-season Anopheles coluzzii in the Sahel. Evol. Appl. 10, 704–717 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
18.Huestis, D. L. et al. Windborne long-distance migration of malaria mosquitoes in the Sahel. Nature 574, 404–408 (2019).CAS 
PubMed 

Google Scholar 
19.Wondji, C., Simard, F. & Fontenille, D. Evidence for genetic differentiation between the molecular forms M and S within the Forest chromosomal form of Anopheles gambiae in an area of sympatry. Insect Mol. Biol. 11, 11–19 (2002).
CAS 
PubMed 

Google Scholar 
20.Simard, F., Nchoutpouen, E., Toto, J. C. & Fontenille, D. Geographic distribution and breeding site preference of Aedes albopictus and Aedes aegypti (Diptera: Culicidae) in Cameroon, Central Africa. J. Med. Entomol. 42, 726–731 (2005).PubMed 

Google Scholar 
21.Roux, O., Diabaté, A. & Simard, F. Divergence in threat sensitivity among aquatic larvae of cryptic mosquito species. J. Anim. Ecol. 83, 702–711 (2014).PubMed 

Google Scholar 
22.Costantini, C. et al. Living at the edge: Biogeographic patterns of habitat segregation conform to speciation by niche expansion in Anopheles gambiae. BMC Ecol. 9 (2009).23.The Anopheles gambiae 1000 Genomes Consortium. Genetic diversity of the African malaria vector Anopheles gambiae. Nature 552, 96–100 (2017).PubMed Central 

Google Scholar 
24.Oliva, C. F., Benedict, M. Q., Lempérière, G. & Gilles, J. Laboratory selection for an accelerated mosquito sexual development rate. Malar. J. 10, 135 (2011).PubMed 
PubMed Central 

Google Scholar 
25.Munhenga, G. et al. Evaluating the potential of the sterile insect technique for malaria control: Relative fitness and mating compatibility between laboratory colonized and a wild population of Anopheles arabiensis from the Kruger National Park, South Africa. Parasit. Vectors 4, 208 (2011).PubMed 
PubMed Central 

Google Scholar 
26.Lee, H. L. et al. Mating compatibility and competitiveness of transgenic and wild type Aedes aegypti (L.) under contained semi-field conditions. Transgenic Res. 22, 47–57 (2013).CAS 
PubMed 

Google Scholar 
27.Damiens, D. et al. Cross-Mating compatibility and competitiveness among Aedes albopictus strains from distinct geographic origins-implications for future application of sit programs in the south west Indian ocean islands. PLoS One 11, e0163788 (2016).PubMed 
PubMed Central 

Google Scholar 
28.Zheng, X. et al. Incompatible and sterile insect techniques combined eliminate mosquitoes. Nature 572, 56–61 (2019).CAS 
ADS 

Google Scholar 
29.Aguilar, R. et al. Genome-wide analysis of transcriptomic divergence between laboratory colony and field Anopheles gambiae mosquitoes of the M and S molecular forms. Insect Mol. Biol. 19, 695–705 (2010).CAS 
PubMed 
PubMed Central 

Google Scholar 
30.Sawadogo, P. S. et al. Swarming behaviour in natural populations of Anopheles gambiae and An. coluzzii: Review of 4 years survey in rural areas of sympatry, Burkina Faso (West Africa). Acta Trop. 130, 24–34 (2014).
Google Scholar 
31.Poda, S. B. et al. Sex aggregation and species segregation cues in swarming mosquitoes: Role of ground visual markers. Parasit. Vectors 12, 589 (2019).PubMed 
PubMed Central 

Google Scholar 
32.Ekechukwu, N. E. et al. Heterosis increases fertility, fecundity, and survival of laboratory-produced F1 hybrid males of the malaria mosquito Anopheles coluzzii. G3 Genes Genomes Genet. 5, 2693–2709 (2015).CAS 

Google Scholar 
33.Ng’habi, K. R. et al. Colonization of malaria vectors under semi-field conditions as a strategy for maintaining genetic and phenotypic similarity with wild populations. Malar. J. 14, 10 (2015).PubMed 
PubMed Central 

Google Scholar 
34.Huho, B. J. et al. Nature beats nurture: A case study of the physiological fitness of free-living and laboratory-reared male Anopheles gambiae s.l. J. Exp. Biol. 210, 2939–2947 (2007).CAS 
PubMed 

Google Scholar 
35.Ferguson, H. M., John, B., Ng, K. & Knols, B. G. J. Redressing the sex imbalance in knowledge of vector biology. Trends Ecol. Evol. 20, 202–209 (2005).PubMed 

Google Scholar 
36.Hassan, M., El-Motasim, W. M., Ahmed, R. T. & El-Sayed, B. B. Prolonged colonisation, irradiation, and transportation do not impede mating vigour and competitiveness of male Anopheles arabiensis mosquitoes under semi-field conditions in Northern Sudan. Malar. World J. 1 (2010).37.Yamada, H., Vreysen, M. J. B., Gilles, J. R. L., Munhenga, G. & Damiens, D. D. The effects of genetic manipulation, dieldrin treatment and irradiation on the mating competitiveness of male Anopheles arabiensis in field cages. Malar. J. 13, 1–10 (2014).
Google Scholar 
38.Munhenga, G. et al. Mating competitiveness of sterile genetic sexing strain males (GAMA) under laboratory and semi-field conditions : Steps towards the use of the Sterile Insect Technique to control the major malaria vector Anopheles arabiensis in South Africa. Parasit. Vectors 9, 1–12 (2016).
Google Scholar 
39.Assogba, B. S. et al. Characterization of swarming and mating behaviour between Anopheles coluzzii and Anopheles melas in a sympatry area of Benin. Acta Trop. 132S, 1–11 (2013).
Google Scholar 
40.Charlwood, J. D. et al. The swarming and mating behaviour of Anopheles gambiae s.s. (Diptera: Culicidae) from São Tomé Island. J. Vector Ecol. 27, 178–183 (2002).CAS 
PubMed 

Google Scholar 
41.Diabate, A. et al. Natural swarming behaviour of the molecular M form of Anopheles gambiae. Trans. R. Soc. Trop. Med. Hyg. 97, 713–716 (2003).CAS 
PubMed 

Google Scholar 
42.Manoukis, N. C. et al. Structure and dynamics of male swarms of Anopheles gambiae. J. Med. Entomol. 46, 227–235 (2009).PubMed 

Google Scholar 
43.Aldersley, A. et al. Too ‘sexy’ for the field? Paired measures of laboratory and semi-field performance highlight variability in the apparent mating fitness of Aedes aegypti transgenic strains. Parasit. Vectors 12, 357 (2019).PubMed 
PubMed Central 

Google Scholar 
44.Pantoja-Sánchez, H., Gomez, S., Velez, V., Avila, F. W. & Alfonso-Parra, C. Precopulatory acoustic interactions of the New World malaria vector Anopheles albimanus (Diptera: Culicidae). Parasit. Vectors 12, 386 (2019).PubMed 
PubMed Central 

Google Scholar 
45.Gibson, G., Warren, B. & Russell, I. J. Humming in tune: Sex and species recognition by mosquitoes on the wing. JARO 540, 527–540 (2010).
Google Scholar 
46.Pennetier, C., Warren, B., Dabiré, K. R., Russell, I. J. & Gibson, G. ‘Singing on the wing’ as a mechanism for species recognition in the malarial mosquito Anopheles gambiae. Curr. Biol. 20, 131–136 (2010).CAS 
PubMed 

Google Scholar 
47.Caputo, B. et al. Comparative analysis of epicuticular lipid profiles of sympatric and allopatric field populations of Anopheles gambiae s.s. molecular forms and An. arabiensis from Burkina Faso (West Africa). Insect Biochem. Mol. Biol. 37, 389–398 (2007).CAS 
PubMed 

Google Scholar 
48.Ferguson, H. M. & Read, A. F. Genetic and environmental determinants of malaria parasite virulence in mosquitoes. Proc. R. Soc. B Biol. Sci. 269, 1217–1224 (2002).CAS 

Google Scholar 
49.Niang, A. et al. Semi-field and indoor setups to study malaria mosquito swarming behavior. Parasit. Vectors 12, 446 (2019).PubMed 
PubMed Central 

Google Scholar 
50.Santolamazza, F. et al. Insertion polymorphisms of SINE200 retrotransposons within speciation islands of Anopheles gambiae molecular forms. Malar. J. 7, 163 (2008).
PubMed 
PubMed Central 

Google Scholar 
51.Vantaux, A. et al. Larval nutritional stress affects vector life history traits and human malaria transmission. Sci. Rep. 6, 36778 (2016).CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
52.Crawley, M. J. The R Book. (Ltd, Sons, 2012). https://doi.org/10.1002/9780470515075.53.Hothorn, T., Bretz, F., Westfall, P. & Heiberger, R. M. Package ‘multcomp’ title simultaneous inference in general parametric models. Biom. J. 50, 346–363 (2016).
Google Scholar  More