Ecology

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

General characterization of seven newly isolated HMO-2011-type phagesIn this study, we used four Roseobacter strains (FZCC0040, FZCC0042, FZCC0012, and FZCC0089) and one SAR11 strain (HTCC1062) to isolate phages. FZCC0040 and FZCC0042 belong to the Roseobacter RCA lineage [22], FZCC0012 shares 99.8% 16S rRNA gene identity with Roseobacter strain HIMB11 [57], and FZCC0089 belongs to a newly identified Roseobacter lineage located close to HIMB11 and SAG-019 lineages (Supplementary Fig. 1).A total of seven phages were newly isolated and analyzed in this study (Table 1). The complete phage genomes range in size from 52.7 to 54.9 kb, harbor 62 to 84 open reading frames (ORFs), and feature a G + C content ranging from 33.8 to 48.6%. Compared to other HMO-2011-type phages, pelagiphage HTVC033P has a relatively lower G + C content of 33.8%, similar to the G + C content of its host HTCC1062 (29.0%) and of other described pelagiphages [21, 26,27,28]. The G + C content of other six roseophages ranges from 42.2 to 48.6%, which is also similar to the G + C content of the hosts they infect (44.8 to 54.1%).Despite their distinct host origins, these phage genomes show considerable similarity in terms of gene content and genome architecture (Fig. 1). They all display clear similarity with the previously reported SAR116 phage HMO-2011 [20] and HMO-2011-type RCA phages [22]. Overall, these phages share 19.2 to 79.1% of their genes with previously reported HMO-2011-type phages and all contain homologues of HMO-2011-type DNA replication and metabolism genes, structural genes, and DNA packaging genes. Moreover, their overall genome structure is conserved with that of HMO-2011-type phages. Considering these observations, we tentatively classified these seven phages into the HMO-2011-type group. Of the 11 currently known HMO-2011-type isolates, one infects the SAR116 strain IMCC1322, one infects the SAR11 strain HTCC1062, and the remaining nine all infect Roseobacter strains; this suggest that HMO-2011-type phages infect diverse bacterial hosts. HTVC033P is the first pelagiphage identified to belong to the HMO-2011-type viral group. Our study has also increased the number of known types of pelagiphages. To date, pelagiphages belonging to a total of nine distinct viral groups have been isolated and analyzed [21, 26,27,28].Fig. 1: Alignment and comparison of genomes of HMO-2011-type isolates and representative HMO-2011-type MVGs from major subgroups.HMO-2011-type phage isolates are shown in red. Phages isolated in this study are indicated with red asterisks. Predicted open reading frames (ORFs) are represented by arrows, with the left or right arrow points indicating the direction of their transcription. The numbers inside the arrows indicate ORF numbers. ORFs annotated with known functions are marked using distinct colors according to their functions. HMO-2011-type core genes are indicated with blue asterisks. The color of the shading connecting homologous genes indicates the level of amino acid identity between the genes. To clearly present the genomic comparison, several MVGs were rearranged to start from the same gene as in the HMO-2011-type phages. DNAP DNA polymerase, Endo endonuclease, RNR ribonucleoside-triphosphate reductase, PhoH phosphate starvation-inducible protein, MazG MazG nucleotide pyrophosphohydrolase domain protein, ThyX thymidylate synthase, GRX glutaredoxin, TerS terminase small subunit, TerL terminase large subunit.Full size imageIdentification and sequence analyses of HMO-2011-type MVGsTo identify HMO-2011-type MVGs, we performed a metagenomic mining and retrieved a total of 207 HMO-2011-type MVGs (≥50% genome completeness) from viromes in the worldwide ocean, from tropical to polar oceans (Supplementary Table 1). These MVGs range in size from 29.2 to 67.9 kb and their G + C content range from 31.3 to 52.4%. In addition, 45 HMO-2011-type MVGs were also identified from some non-marine habitats, suggesting that HMO-2011-type phages are widely distributed worldwide (Supplementary Table 1).Genomic analysis confirmed that all HMO-2011-type MVGs exhibit genomic synteny with HMO-2011-type phages (Fig. 1). Although some of these HMO-2011-type MVGs are highly similar to their cultivated relatives, most MVGs appear to have more genomic variations. To resolve the evolutionary relationship among the HMO-2011-type phages, a phylogenetic tree was constructed based on the concatenated sequences of five core genes. We found that HMO-2011-type phages are evolutionarily diverse and can be separated into at least 10 well-supported subgroups ( >2 members), with 140 MVGs clustering into previously identified HMO-2011-type groups (subgroups I and III in Fig. 2A) [22], and the remaining 67 MVGs forming new subgroups (Fig. 2A). Among these HMO-2011-type subgroups, three contain cultivated representatives (subgroups I, III, and IX). Subgroup I contains the greatest number of phages, including six cultivated representatives and 123 MVGs (Fig. 2A). The cultivated representatives in subgroup I include a phage that infect SAR116 strain and five phages that infect Roseobacter strains. Subgroup III contains four cultivated representatives that infect two Roseobacter strains, and 17 MVGs. Pelagiphage HTVC033P and nine MVGs form subgroup IX. Other subgroups have no cultivated representatives yet. The results of phylogenomic analysis showed that subgroups I to VI are closely related, whereas subgroups VII to X are located on a separate branch and are more distinct from the subgroups I to VI, which suggests that these subgroups are more evolutionarily distant. A phylogenomic-based approach with GL-UVAB workflow [53] was also performed to cluster these HMO-2011-type genomes, which showed similar grouping results (Supplementary Fig. 2).Fig. 2: Phylogenomic and shared-gene analyses of HMO-2011-type phages.A A maximum-likelihood tree was constructed using concatenated sequences of five hallmark genes. HMO-2011-type phages were grouped into 10 subgroups based on the phylogeny. Shading is used to indicate the subgroups. HMO-2011-type phage isolates are shown in red. Genomes containing an integrase gene are indicated by red triangles. The G + C content and completeness of the genomes are indicated. Scale bar indicates the number of amino acid substitutions per site. B Heatmap showing the percentage of shared genes between HMO-2011-type genomes. Phages in the same subgroup are boxed.Full size imageA previous study suggested the use of the percentage of shared proteins as a means of defining phage taxonomic ranks and proposed that phages with ≥20 and ≥40% orthologous proteins in common can be grouped at the taxonomic ranks of subfamily and genus, respectively [58]. Overall, most of the calculated percentages between HMO-2011-type genomes fall within the 20 to 100% range and most of the percentages between genomes within the same subgroup fall within the 40 to 100% range (Fig. 2B). Therefore, our results suggest that the HMO-2011-type is roughly a subfamily-level phage taxonomic group containing at least ten genus-level subgroups in the Podoviridae family.Conserved genomic structure and variation in HMO-2011-type phagesOf the 1235 orthologous protein groups (≥2 members) identified in HMO-2011-type genomes, only 254 proteins groups could be assigned putative biological functions (Supplementary Table 2). Comparative genomic analysis clearly revealed the conserved functional module structure of all HMO-2011-type genomes. All HMO-2011-type phage genomes can be roughly divided into the DNA metabolism and replication module, structural module and DNA packaging module (Fig. 1). Most of the homologous genes are scattered in similar loci of the HMO-2011-type genomes. Core genome analysis based on complete HMO-2011-type genomes revealed that HMO-2011-type genomes share a common set of ten core genes (Fig. 1). These core genes are mostly genes related to essential function in phage replication and development, including genes encoding DNA helicase, DNA primase, DNA polymerase (DNAP), portal protein, capsid protein, and terminase small and large subunits (TerL and TerS) as well as several genes with no known function, suggesting that phages in this group employ similar overall infection and propagation processes (Fig. 1).Most members in subgroups I and III and one member in subgroup II possess a tyrosine integrase gene (int) located upstream of the DNA replication and metabolism module, whereas all subgroup IV to X genomes contain no identifiable lysogeny-related genes. This result suggests that members of subgroups IV to X might be obligate lytic phages. Integrase genes typically occur in the genomes of temperate phages and are responsible for site-specific recombination between phage and host bacterial genomes [59, 60]. In subgroup III, RCA phage CRP-3 has been experimentally demonstrated to be capable of integrating into the host genome [22]. Thus, certain int-containing HMO-2011-type phages are also likely to be temperate phages.In the DNA metabolism and replication modules, genes encoding DNA primase, DNA helicase, DNAP, ribonucleotide reductase (RNR), and endonuclease can be identified; and DNA helicase, DNA primase, and DNAP are core to all HMO-2011-type phages. All reported HMO-2011-type phages contain an atypical DNAP, in which a partial DnaJ central domain is located between the exonuclease domain and the DNA polymerase domain [20, 22]. The Escherichia coli DnaJ protein, a co-chaperone [61], has been shown to be involved in diverse functions [62] and to be critical for the replication of phage Lambda [63,64,65]. The sequence analysis revealed that DNAP sequences of these seven new HMO-2011-type phages and 207 MVGs also present this unusual domain structure and contain two repeats of the CXXCXGXG motifs involved in zinc binding [66] in the partial DnaJ domain (Supplementary Fig. 3). RNR gene is frequently detected in subgroups I, II, III, IV, V, and X genomes but not in the other subgroup genomes. RNRs, which are widely distributed in diverse phage genomes, are involved in catalyzing the reduction of ribonucleotides to deoxyribonucleotides, and thus play a crucial role in providing deoxyribonucleoside triphosphates for phage DNA biosynthesis and repair [67,68,69]. RNR genes clustered with the RNR gene in phage HMO-2011 were previously reported to dominate the class II viral RNRs in examined marine viromes [69]. In the remaining two modules, genes involved in phage structure (e.g., genes encoding capsid and portal proteins), packaging of DNA (TerL and TerS genes), and cell lysis were detected. The proteins encoded by these genes play key roles in phage morphogenesis and virion release.Examination of the distribution of the orthologous groups among the subgroups revealed clear pan-genome differences in various subgroups (Fig. 3). Most subgroups harbor subgroup-specific genes not identified in other subgroups, although  no function has yet been assigned to most of these genes. Notably, the phages in subgroups VII, VIII, and IX possess genomic features that differentiate them from phages in other subgroups, specifically with regard to the G + C content and gene content. The members of these three subgroups are closely related to each other in the phylogenetic tree and harbor several subgroup-specific genes. The G + C content of the phage genomes in these subgroups ranges from 31.9 to 35.4%, significantly smaller than other subgroups but similar to the G + C content of SAR11 bacteria and other known pelagiphages. HTVC033P is the only cultivated representative of subgroup IX. The aforementioned results suggest that the phages in subgroup VII, VIII, and IX might have related bacterial hosts and are highly likely to be pelagiphages. The host prediction using RaFAH tool also assigned Pelagibacter as their potential hosts (Supplementary Table 1). Subgroup X is located near these three subgroups in the phylogenetic tree, and the G + C content of the phages in this subgroup ranges from 34.4 to 39.0%. The host prediction assigned Roseobacter as their potential hosts. The hosts of this subgroup still remain to be experimentally investigated.Fig. 3: Distribution and functional classification of orthologous protein groups across HMO-2011-type genomes.Only orthogroups containing >10 members or showing subgroup-specific features are shown. Subgroup-specific genes are boxed in red. Genes that are absent in a specific subgroup are boxed in orange.Full size imageMetabolic capabilities of HMO-2011-type phagesAll HMO-2011-type phage genomes harbor several host-derived auxiliary metabolic genes (AMGs) potentially involved in diverse metabolic processes. Some AMGs in HMO-2011-type phages have been discussed previously [20, 22].Subgroups VII, VIII, IX, and X possess distinct AMGs as compared with the other subgroups. For example, the genes encoding FAD-dependent thymidylate synthase (ThyX, PF02511) and MazG pyrophosphohydrolase domains are absent in all subgroups VII, VIII, IX, and X genomes but frequently detected in other subgroup genomes. ThyX protein is essential for the conversion of dUMP to dTMP mediated by an FAD coenzyme and is therefore a key enzyme involved in DNA synthesis [70, 71]. The thyX gene is commonly found in microbial genomes and phage genomes. Phage-encoded ThyX has been suggested to compensate for the loss of host-encoded ThyA and thus play crucial roles in phage nucleic acid synthesis and metabolism during infection [72]. Except in the case of subgroups VII, VIII, IX, and X genomes, the mazG gene, which encodes a nucleoside triphosphate pyrophosphohydrolase is sporadically distributed in HMO-2011-type genomes. MazG protein is predicted to be a regulator of nutrient stress and programmed cell death [73] and has been hypothesized to promote phage survival by keeping the host alive during phage propagation [74]. The Escherichia coli MazG can interfere with the function of the MazEF toxin–antitoxin system by decreasing the cellular level of (p)ppGpp [73]. However, a recent study showed that a cyanophage MazG has no binding or hydrolysis activity against alarmone (p)ppGpp but has high hydrolytic activity toward dGTP and dCTP, and it was speculated to play a role in hydrolyzing high G + C host genome for phage replication [75]. Whether the MazG proteins encoded by HMO-2011-type phages play a similar role in phage propagation remained to be investigated.Five MVGs in subgroup I contain a gene encoding a DraG-like family ADP-ribosyl hydrolase (ARH). In cellular ADP-ribosylation systems, ARH catalyzes the cleavage of the ADP-ribose moiety, and thereby counteract the effects of ADP-ribosyl transferases [76]. It has been reported that ARH in Rhodospirillum rubrum regulates the nitrogen fixation [77]. However, the function of this phage-encoded ARH in the phage propagation process remains unclear.We also observed that several MVGs possess genes involved in iron–sulfur (Fe–S) cluster biosynthesis, including an Fe–S cluster assembly scaffold gene (iscU) that involved in Fe–S cluster assembly and transfer [78] and an Fe–S cluster insertion protein gene (erpA). Fe–S cluster participates in a wide variety of cellular biological processes [79]. The discovery of these genes suggests that these phages may play important roles in Fe–S cluster biogenesis and function.The gene encoding sodium-dependent phosphate transport protein (PF02690) has been identified in eight subgroup I genomes. The Na/Pi cotransporter family protein is responsible for high-affinity, sodium-dependent Pi uptake, and thus the protein plays a critical role in maintaining phosphate homeostasis [80]. This gene might function in the transport of phosphate into cells during phage infection. The presence of Na/Pi cotransporter genes suggests that some HMO-2011-type phages may have the potential to regulate host phosphate uptake in phosphate-limited ocean environments in order to benefit phage replication and propagation.Identification and phylogenetic analysis of HMO-2011-type DNAPsThe genetic diversity and geographically distribution of HMO-2011-type phages in marine environments was further inferred from DNAP gene analyses. A total of 2433 HMO-2011-type DNAP sequences with sequence sizes ranging from 540 to 779 amino acids were identified and subjected to phylogenetic analysis (Supplementary Table 3).Among the identified HMO-2011-type DNAPs, 2030 sequences were retrieved from the GOV 2.0 Tara expedition upper-ocean viral populations (0–1000 m), from tropical to polar regions. HMO-2011-type DNAP genes were identified from all analyzed upper-ocean viromes, suggesting the global prevalence of HMO-2011-type phages in upper oceans.A previous study revealed that marine viromes contain various types of tailed phage genomes that encode a family A DNAP gene [81]. To estimate the importance of HMO-2011-type phages, we calculated the proportion of HMO-2011-type DNAPs based on the number of HMO-2011-type DNAP sequences and the total number of family A DNAP sequences ( >470 aa) in each GOV 2.0 viral population dataset. This analysis revealed that HMO-2011-type DNAPs accounted for up to 19.7% of all family A DNAPs in each GOV 2.0 dataset (Supplementary Table 4). We found that the HMO-2011-type DNAP sequences appear to be more dominant in epipelagic viromes than in mesopelagic viromes (p  More

1.Caddy, J. F. & Gulland, J. A. Historical patterns of fish stocks. Mar. Policy 7, 267–278 (1983).
Google Scholar 
2.Steele, J. H. & Henderson, E. W. Coupling between physical and biological scales. Philos. Trans. R. Soc. Lond. Ser. B-Biol. Sci. 343, 5–9 (1994).
Google Scholar 
3.Bjornstad, O. N., Fromentin, J. M., Stenseth, N. C. & Gjosaeter, J. Cycles and trends in cod populations. Proc. Natl Acad. Sci. USA 96, 5066–5071 (1999).CAS 
PubMed 
PubMed Central 

Google Scholar 
4.Piatt, J. F. et al. Extreme mortality and reproductive failure of common murres resulting from the northeast Pacific marine heatwave of 2014–2016. PLoS One 15, e0226087 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
5.Oremus, K. L. Climate variability reduces employment in New England fisheries. Proc. Natl Acad. Sci. USA 16, 26444–26449 (2018).
Google Scholar 
6.Shelton, A. O. & Mangel, M. Fluctuations of fish populations and the magnifying effect of fishing. Proc. Natl Acad. Sci. USA 108, 7075–7080 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
7.Essington, T. E. et al. Fishing amplifies forage fish population collapses. Proc. Natl Acad. Sci. USA 112, 6648–6652 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
8.Memarzadeha, M., Britten, G. L., Wormd, B. & Boettigere, C. Rebuilding global fisheries under uncertainty. Proc. Natl Acad. Sci. USA 116, 15985–15990 (2019).
Google Scholar 
9.Pauly, D. & Zeller, D. Sea Around Us Concepts, Design and Data (seaaroundus.org) (2015).10.Bjornstad, O. N., Nisbet, R. M. & Fromentin, J. M. Trends and cohort resonant effects in age-structured populations. J. Anim. Ecol. 73, 1157–1167 (2004).
Google Scholar 
11.Botsford, L. W., Holland, M. D., Field, J. C. & Hastings, A. Cohort resonance: a significant component of fluctuations in recruitment, egg production, and catch of fished populations. ICES J. Mar. Sci. 71, 2158–2170 (2014).
Google Scholar 
12.Di Lorenzo, E. & Ohman, M. D. A double-integration hypothesis to explain ocean ecosystem response to climate forcing. Proc. Natl Acad. Sci. USA 110, 2496–2499 (2013).PubMed 
PubMed Central 

Google Scholar 
13.Bjorkvoll, E. et al. Stochastic population dynamics and life-history variation in marine fish species. Am. Naturalist 180, 372–387 (2012).
Google Scholar 
14.Hsieh, C. H. et al. Fishing elevates variability in the abundance of exploited species. Nature 443, 859–862 (2006).CAS 
PubMed 

Google Scholar 
15.Beamish, R. J., McFarlane, G. A. & Benson, A. Longevity overfishing. Prog. Oceanogr. 68, 289–302 (2006).
Google Scholar 
16.Anderson, C. N. K. et al. Why fishing magnifies fluctuations in fish abundance. Nature 452, 835–839 (2008).CAS 
PubMed 

Google Scholar 
17.Hutchings, J. A. & Myers, R. A. Effect of age on the seasonality of maturation and spawning of Atlantic cod, Gadus morhua, in the northwest Atlantic. Can. J. Fish. Aquat. Sci. 50, 2468–2474 (1993).
Google Scholar 
18.Bobko, S. J. & Berkeley, S. A. Maturity, ovarian cycle, fecundity, and age-specific parturition of black rockfish (Sebastes melanops). Fish. Bull. 102, 418–429 (2004).
Google Scholar 
19.Berkeley, S. A., Chapman, C. & Sogard, S. M. Maternal age as a determinant of larval growth and survival in a marine fish, Sebastes melanops. Ecology 85, 1258–1264 (2004).
Google Scholar 
20.Longhurst, A. Murphy’s law revisited: longevity as a factor in recruitment to fish populations. Fish. Res. 56, 125–131 (2002).
Google Scholar 
21.Stawitz, C. C. & Essington, T. E. Somatic growth contributes to population variation in marine fishes. J. Anim. Ecol. 88, 315–329 (2019).PubMed 

Google Scholar 
22.Estes, J. A. et al. Trophic downgrading of planet Earth. Science 333, 301–306 (2011).CAS 
PubMed 

Google Scholar 
23.Hollowed, A. B., Hare, S. R. & Wooster, W. S. Pacific basin climate variability and patterns of Northeast Pacific marine fish production. Prog. Oceanogr. 49, 257–282 (2001).
Google Scholar 
24.Holsman, K. K., Aydin, K., Sullivan, J., Hurst, T. & Kruse, G. H. Climate effects and bottom-up controls on growth and size-at-age of Pacific halibut (Hippoglossus stenolepis) in Alaska (USA). Fish. Oceanogr. 28, 345–358 (2019).
Google Scholar 
25.Whitten, A. R., Klaer, N. L., Tuck, G. N. & Day, R. W. Accounting for cohort-specific variable growth in fisheries stock assessments: A case study from south-eastern Australia. Fish. Res. 142, 27–36 (2013).
Google Scholar 
26.Heessen, H. J. L., Daan, N. & Ellis, J. R. Fish atlas of the Cebtic Sea, North Sea, and Baltic Sea (KNNV Publishing and Wageningen Academic Publishers, 2015).27.Froese, R. & Pauly, D. FishBase, version (01/2021) https://www.fishbase.org (2021).28.Munch, S. B. & Salinas, S. Latitudinal variation in lifespan within species is explained by the metabolic theory of ecology. Proc. Natl Acad. Sci. USA 106, 13860–13864 (2009).CAS 
PubMed 
PubMed Central 

Google Scholar 
29.Beukhof, E. et al. Marine fish traits follow fast-slow continuum across oceans. Sci. Rep. 9, 17878 (2019).30.Pauly, D. On the interrelationships between natural mortality, growth parameters, and mean environmental temperature in 175 fish stocks. J. Cons. int. Explor. Mer. 39, 175–192 (1980).
Google Scholar 
31.Audzijonyte, A. et al. Fish body sizes change with temperature but not all species shrink with warming. Nat. Ecol. Evol. 4, 1–6 (2020).
Google Scholar 
32.Audzijonyte, A. et al. Is oxygen limitation in warming waters a valid mechanism to explain decreased body sizes in aquatic ectotherms? Glob. Ecol. Biogeogr. 28, 64–77 (2019).
Google Scholar 
33.Forster, J., Hirst, A. G. & Atkinson, D. Warming-induced reductions in body size are greater in aquatic than terrestrial species. Proc. Natl Acad. Sci. USA 109, 19310–19314 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
34.Block, B. A. et al. Tracking apex marine predator movements in a dynamic ocean. Nature 475, 86–90 (2011).CAS 
PubMed 

Google Scholar 
35.Behrenfeld, M. J. & Falkowski, P. G. Photosynthetic rates derived from satellite-based chlorophyll concentration. Limnol. Oceanogr. 42, 1–20 (1997).CAS 

Google Scholar 
36.Ives, A. R. Measuring resilience in stochastic-systems. Ecol. Monogr. 65, 217–233 (1995).
Google Scholar 
37.Alheit, J. & Niquen, M. Regime shifts in the Humboldt Current ecosystem. Prog. Oceanogr. 60, 201–222 (2004).
Google Scholar 
38.Pinsky, M. L., Jensen, O. P., Ricard, D. & Palumbi, S. R. Unexpected patterns of fisheries collapse in the world’s oceans. Proc. Natl Acad. Sci. USA 108, 8317–8322 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
39.Spencer, P. D. & Collie, J. S. Patterns of population variability in marine fish stocks. Fish. Oceanogr. 6, 188–204 (1997).
Google Scholar 
40.FAO. The State of World Fisheries and Aquaculture 2018 – Meeting the Sustainable Development Goals. (Food and Agriculture Organization of the United Nations, Rome, 2018).41.Barnett, L. A. K., Branch, T. A., Ranasinghe, R. A. & Essington, T. E. Old-growth fishes become scarce under fishing. Curr. Biol. 27, 2843–2848 (2017).CAS 
PubMed 

Google Scholar 
42.Rouyer, T. et al. Shifting dynamic forces in fish stock fluctuations triggered by age truncation? Glob. Change Biol. 17, 3046–3057 (2011).
Google Scholar 
43.Coumou, D. & Rahmstorf, S. A decade of weather extremes. Nat. Clim. Change 2, 491–496 (2012).
Google Scholar 
44.Easterling, D. R. et al. Climate extremes: Observations, modeling, and impacts. Science 289, 2068–2074 (2000).CAS 
PubMed 

Google Scholar 
45.Portner, H. O. & Peck, M. A. Climate change effects on fishes and fisheries: towards a cause-and-effect understanding. J. Fish. Biol. 77, 1745–1779 (2010).CAS 
PubMed 

Google Scholar 
46.Pinsky, M. L., Worm, B., Fogarty, M. J., Sarmiento, J. L. & Levin, S. A. Marine taxa track local climate velocities. Science 341, 1239–1242 (2013).CAS 
PubMed 

Google Scholar 
47.de Gee, A. & Kikkert, A. H. Analysis of the grey gurnard (Eutrigla gurnardus) samples collected during the 1991 international stomach sample project. ICES Document CM 1993/G:14, 25 (1993).48.Sparholt, H. In Fish Atlas of the Celtic Sea, North Sea, and Baltic Sea (eds Heessen, H., Daan, N., & Ellis, J. R.) 377–381 (KNNV Publishiing and Wageningen Academic Publishers, 2015).49.Arnott, S. A. & Ruxton, G. D. Sandeel recruitment in the North Sea: demographic, climate and trophic effects. Mar. Ecol. Prog. Ser. 238, 199–210 (2002).
Google Scholar 
50.van Deurs, M., van Hal, R., Tomczak, M. T., Jonasdottir, S. H. & Dolmer, P. Recruitment of lesser sandeel Ammodytes marinus in relation to density dependence and zooplakton composition. Mar. Ecol. Prog. Ser. 381, 249–258 (2009).
Google Scholar 
51.Capuzzo, E. et al. A decline in primary production in the North Sea over 25 years, associated with reductions in zooplankton abundance and fish stock recruitment. Glob. Change Biol. 24, E352–E364 (2018).
Google Scholar 
52.Rayner, N. A. et al. Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. J. Geophys. Res. D: Atmospheres 108, ACL 2-1–ACL 2–29 (2003).
Google Scholar 
53.Papworth, D. J., Marini, S. & Conversi, A. Novel, unbiased analysis approach for investigating population dynamics: A case study on Calanus finmarchicus and its decline in the North Sea. PLoS One 11, e0158230 (2016).PubMed 
PubMed Central 

Google Scholar 
54.Bergstad, O. A., Hoines, A. S. & Jorgensen, T. Growth of sandeel Ammodytes marinus, in the northern North Sea and Norwegian coastal waters. Fish. Res. 56, 9–23 (2002).
Google Scholar 
55.Wright, P. J. Otolith microstructure of the lesser sandeel, Ammodytes marinus. J. Mar. Biol. Assoc. U.K. 73, 245–248 (1993).
Google Scholar 
56.Sell, A. & Heessen, H. in Fish atlas of the Celtic Sea, North Sea, and Baltic Sea (eds Heessen, H., Daan, N., & Ellis, J. R.) 295−299 (KNNV Publishing and Wageningen Academic Publishers, 2015).57.Bergstad, O. A., Hoines, A. S. & Kruger-Johnsen, E. M. Spawning time, age and size at maturity, and fecundity of sandeel, Ammodytes marinus, in the north-eastern North Sea and in unfished coastal waters off Norway. Aquat. Living Resour. 14, 293–301 (2001).
Google Scholar 
58.Pyper, B. J. & Peterman, R. M. Comparison of methods to account for autocorrelation in correlation analyses of fish data. Can. J. Fish. Aquat. Sci. 55, 2127–2140 (1998).
Google Scholar 
59.van der Sleen, P. et al. Non-stationary responses in anchovy (Engraulis encrasicolus) recruitment to coastal upwelling in the Southern Benguela. Mar. Ecol. Prog. Ser. 596, 155–164 (2018).
Google Scholar 
60.Cushing, D. H. Upwelling and production on fish. Adv. Mar. Biol. 9, 255–334 (1971).
Google Scholar 
61.Pauly, D. & Lam, V. W. Y. In Large marine ecosystems: Status and Trends (eds IOC-UNESCO and UNEP) 113–137 (United Nations Environmental Programme, 2016). More

1.Rodrigues, A. S. B. et al. New mitochondrial and nuclear evidences support recent demographic expansion and an atypical phylogeographic pattern in the spittlebug Philaenus spumarius (Hemiptera, Aphrophoridae). PLoS ONE 9, 1–12. https://doi.org/10.1371/journal.pone.0098375 (2014).CAS 
Article 

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

Google Scholar 
3.Saponari, M. et al. Infectivity and transmission of Xylella fastidiosa by Philaenus spumarius (Hemiptera: Aphrophoridae) in Apulia Italy. J. Econ. Entomol. 107, 1316–1319. https://doi.org/10.1603/EC14142 (2014).Article 
PubMed 

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

Google Scholar 
5.Cavalieri, V. et al. Transmission of Xylella fastidiosa subspecies pauca sequence type 53 by different insect species. Insects 10, 324. https://doi.org/10.3390/insects10100324 (2019).Article 
PubMed Central 

Google Scholar 
6.Almeida, R. P. P., Blua, M. J., Lopes, J. R. S. & Purcell, A. H. Vector transmission of Xylella fastidiosa: Applying fundamental knowledge to generate disease management strategies. Entomol. Soc. Am. 98, 775–786. https://doi.org/10.1603/0013-8746(2005)098[0775:vtoxfa]2.0.co;2 (2005).Article 

Google Scholar 
7.Schneider, K. et al. Impact of Xylella fastidiosa subspecies pauca in European olives. Proc. Natl. Acad. Sci. U. S. A. 117, 9250–9259. https://doi.org/10.1073/pnas.1912206117 (2020).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
8.Dongiovanni, C. et al. Evaluation of insecticides for the control of juveniles of Philaenus spumarius L., 2015–2017. Arthropod Manag. Tests 43, 2015–2017. https://doi.org/10.1093/amt/tsy073 (2018).Article 

Google Scholar 
9.Fierro, A., Liccardo, A. & Porcelli, F. A lattice model to manage the vector and the infection of the Xylella fastidiosa on olive trees. Sci. Rep. 9, 1–14. https://doi.org/10.1038/s41598-019-44997-4 (2019).CAS 
Article 

Google Scholar 
10.Nyffeler, M. & Benz, G. Spiders in natural pest control: A review. J. Appl. Entomol. 103, 321–339. https://doi.org/10.1111/j.1439-0418.1987.tb00992.x (1987).Article 

Google Scholar 
11Nyffeler, M. & Birkhofer, K. An estimated 400–800 million tons of prey are annually killed by the global spider community. Sci. Nat. https://doi.org/10.1007/s00114-017-1440-1 (2017).Article 

Google Scholar 
12.Nyffeler, M. Ecological impact of spider predation: A critical assessment of Bristowe’s and Turnbull’s estimates. Bull. Br. Arachnol. Soc. 11, 367–373 (2000).
Google Scholar 
13.Phillipson, J. A contribution to the feeding biology of Mitopus morio (F) (Phalangida). J. Anim. Ecol. 29, 35–43. https://doi.org/10.2307/2269 (1960).Article 

Google Scholar 
14Harper, G. & Whittaker, J. The role of natural enemies in the colour polymorphism of Philaenus spumarius (L.). J. Anim. Ecol. 45, 91–104. https://doi.org/10.2307/3769 (1976).Article 

Google Scholar 
15.Benhadi-Marín, J. et al. A guild-based protocol to target potential natural enemies of Philaenus spumarius (Hemiptera: Aphrophoridae), a vector of Xylella fastidiosa (Xanthomonadaceae): A case study with spiders in the olive grove. Insects 11, 100. https://doi.org/10.3390/insects11020100 (2020).Article 
PubMed Central 

Google Scholar 
16.King, R. A., Read, D. S., Traugott, M. & Symondson, W. O. C. Molecular analysis of predation: A review of best practice for DNA-based approaches. Mol. Ecol. 17, 947–963. https://doi.org/10.1111/j.1365-294X.2007.03613.x (2008).CAS 
Article 
PubMed 

Google Scholar 
17.Sint, D., Raso, L., Kaufmann, R. & Traugott, M. Optimizing methods for PCR-based analysis of predation. Mol. Ecol. Resour. 11, 795–801. https://doi.org/10.1111/j.1755-0998.2011.03018.x (2011).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
18.Rejili, M. et al. A PCR-based diagnostic assay for detecting DNA of the olive fruit fly, Bactrocera oleae, in the gut of soil-living arthropods. Bull. Entomol. Res. 106, 695–699. https://doi.org/10.1017/S000748531600050X (2016).CAS 
Article 
PubMed 

Google Scholar 
19Albertini, A. et al. Detection of Bactrocera oleae (Diptera: Tephritidae) DNA in the gut of the soil species Pseudoophonus rufipes (coleoptera: Carabidae). Span. J. Agric. Res. https://doi.org/10.5424/sjar/2018163-12860 (2018).Article 

Google Scholar 
20.Symondson, W. O. C. Molecular identification of prey in predator diets. Mol. Ecol. 11, 627–641. https://doi.org/10.1046/j.1365-294X.2002.01471.x (2002).CAS 
Article 
PubMed 

Google Scholar 
21.Sousa, L. L., Silva, S. M. & Xavier, R. DNA metabarcoding in diet studies: Unveiling ecological aspects in aquatic and terrestrial ecosystems. Environ. DNA 1, 199–214. https://doi.org/10.1002/edn3.27 (2019).Article 

Google Scholar 
22.Juen, A. & Traugott, M. Amplification facilitators and multiplex PCR: Tools to overcome PCR-inhibition in DNA-gut-content analysis of soil-living invertebrates. Soil Biol. Biochem. 38, 1872–1879. https://doi.org/10.1016/j.soilbio.2005.11.034 (2006).CAS 
Article 

Google Scholar 
23.Monzó, C., Sabater-Muñoz, B., Urbaneja, A. & Castańera, P. Tracking medfly predation by the wolf spider, Pardosa cribata Simon, in citrus orchards using PCR-based gut-content analysis. Bull. Entomol. Res. 100, 145–152. https://doi.org/10.1017/S0007485309006920 (2010).CAS 
Article 
PubMed 

Google Scholar 
24.Lantero, E., Matallanas, B., Pascual, S. & Callejas, C. PCR species-specific primers for molecular gut content analysis to determine the contribution of generalist predators to the biological control of the vector of Xylella fastidiosa. Sustainability 10, 4–11. https://doi.org/10.3390/su10072207 (2018).CAS 
Article 

Google Scholar 
25.Cohen, A. C. Extra-oral digestion in predaceous terrestrial Arthropoda. Annu Rev Entomol. 40, 85–103. https://doi.org/10.1146/annurev.en.40.010195.000505 (1995).ADS 
CAS 
Article 

Google Scholar 
26.Krehenwinkel, H., Rödder, N. & Tautz, D. Eco-genomic analysis of the poleward range expansion of the wasp spider Argiope bruennichi shows rapid adaptation and genomic admixture. Glob. Chang. Biol 21, 4320–4332. https://doi.org/10.1111/gcb.13042 (2015).ADS 
Article 
PubMed 

Google Scholar 
27.Kennedy, S. R. et al. High-throughput sequencing for community analysis: the promise of DNA barcoding to uncover diversity, relatedness, abundances and interactions in spider communities. Dev. Genes Evol. 230, 185–201. https://doi.org/10.1007/s00427-020-00652-x (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
28.Hoogendoorn, M. & Heimpel, G. E. PCR-based gut content analysis of insect predators: Using ribosomal ITS-1 fragments from prey to estimate predation frequency. Mol. Ecol. 10, 2059–2067. https://doi.org/10.1046/j.1365-294X.2001.01316.x (2001).CAS 
Article 
PubMed 

Google Scholar 
29.Eitzinger, B., Unger, E. M., Traugott, M. & Scheu, S. Effects of prey quality and predator body size on prey DNA detection success in a centipede predator. Mol. Ecol. 23, 3767–3776. https://doi.org/10.1111/mec.12654 (2014).CAS 
Article 
PubMed 

Google Scholar 
30.Unruh, T. R. et al. Gut content analysis of arthropod predators of codling moth in Washington apple orchards. Biol. Control. 102, 85–92. https://doi.org/10.1016/j.biocontrol.2016.05.014 (2016).Article 

Google Scholar 
31.Rowley, C. et al. PCR-based gut content analysis to identify arthropod predators of Haplodiplosis marginata. Biol. Control 115, 112–118. https://doi.org/10.1016/j.biocontrol.2017.10.003 (2017).CAS 
Article 

Google Scholar 
32.Macías-Hernández, N. et al. Molecular gut content analysis of different spider body parts. PLoS ONE 13, 1–16. https://doi.org/10.1371/journal.pone.0196589 (2018).CAS 
Article 

Google Scholar 
33.Troedsson, C., Simonelli, P., Nägele, V., Nejstgaard, J. C. & Frischer, M. E. Quantification of copepod gut content by differential length amplification quantitative PCR (dla-qPCR). Mar. Biol. 156, 253–259. https://doi.org/10.1007/s00227-008-1079-8 (2009).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
34.Agustí, N., De Vicente, M. C. & Gabarra, R. Development of sequence amplified characterized region (SCAR) markers of Helicoverpa armigera: A new polymerase chain reaction-based technique for predator gut analysis. Mol. Ecol. 8, 1467–1474. https://doi.org/10.1046/j.1365-294X.1999.00717.x (1999).Article 
PubMed 

Google Scholar 
35.Agustí, N., Unruh, T. R. & Welter, S. C. Detecting Cacopsylla pyricola (Hemiptera: Psyllidae) in predator guts using COI mitochondrial markers. Bull. Entomol. Res. 93, 179–185. https://doi.org/10.1079/ber2003236 (2003).Article 
PubMed 

Google Scholar 
36.Aebi, A. et al. Detecting arthropod intraguild predation in the field. Biocontrol 56, 429–440. https://doi.org/10.1007/s10526-011-9378-2 (2011).Article 

Google Scholar 
37.Hosseini, R., Schmidt, O. & Keller, M. A. Factors affecting detectability of prey DNA in the gut contents of invertebrate predators: A polymerase chain reaction-based method. Entomol. Exp. Appl. 126, 194–202. https://doi.org/10.1111/j.1570-7458.2007.00657.x (2008).CAS 
Article 

Google Scholar 
38.Agustí, N. et al. Collembola as alternative prey sustaining spiders in arable ecosystems: Prey detection within predators using molecular markers. Mol Ecol. 12, 3467–3475. https://doi.org/10.1046/j.1365-294X.2003.02014.x (2003).CAS 
Article 
PubMed 

Google Scholar 
39.Greenstone, M. H., Tillman, P. G. & Hu, J. S. Predation of the newly invasive pest Megacopta cribraria (Hemiptera: Plataspidae) in soybean habitats adjacent to cotton by a complex of predators. J Econ Entomol. 107, 947–954. https://doi.org/10.1603/EC13356 (2014).CAS 
Article 
PubMed 

Google Scholar 
40.Welch, K. D., Whitney, T. D. & Harwood, J. D. Non-pest prey do not disrupt aphid predation by a web-building spider. Bull. Entomol. Res. 106, 91–98. https://doi.org/10.1017/S0007485315000875 (2016).CAS 
Article 
PubMed 

Google Scholar 
41.Vincent, J. F. V. Arthropod cuticle: A natural composite shell system. Compos. Part A Appl. Sci. Manuf. 33, 1311–1315. https://doi.org/10.1016/S1359-835X(02)00167-7 (2002).Article 

Google Scholar 
42.Cardoso, P. et al. Rapid biodiversity assessment of spiders (Araneae) using semi-quantitative sampling: A case study in a Mediterranean forest. Insect Conserv. Divers. 1, 71–84. https://doi.org/10.1111/j.1752-4598.2007.00008.x (2008).Article 

Google Scholar 
43.Harwood, J. D., Phillips, S. W., Sunderland, K. D. & Symondson, W. O. C. Secondary predation: Quantification of food chain errors in an aphid-spider-carabid system using monoclonal antibodies. Mol. Ecol. 10, 2049–2057. https://doi.org/10.1046/j.0962-1083.2001.01349.x (2001).CAS 
Article 
PubMed 

Google Scholar 
44.Seabra, S. G. et al. Corrigendum to “Molecular phylogeny and DNA barcoding in the meadow-spittlebug Philaenus spumarius (Hemiptera, Cercopidae) and its related species” [Mol. Phylogenet. Evol. 56 (2010) 462–467]. Mol. Phylogenet. Evol. 152, 106888. https://doi.org/10.1016/j.ympev.2020.106888 (2020).CAS 
Article 
PubMed 

Google Scholar 
45.Folmer, O., Black, M., Hoeh, W., Lutz, R. & Vrijenhoek, R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol. Mar. Biol. Biotechnol. 3, 294–299 (1994).CAS 
PubMed 
PubMed Central 

Google Scholar 
46.Kumar, S., Stecher, G., Li, M., Knyaz, C. & Tamura, K. MEGA X: Molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35, 1547–1549. https://doi.org/10.1093/molbev/msy096 (2018).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
47.Untergasser, A. et al. Primer3-new capabilities and interfaces. Nucleic Acids Res. 40, 1–12. https://doi.org/10.1093/nar/gks596 (2012).CAS 
Article 

Google Scholar 
48.Ye, J. et al. Primer-BLAST: A tool to design target-specific primers for polymerase chain reaction. BMC Bioinform. 13, 134. https://doi.org/10.1186/1471-2105-13-134 (2012).CAS 
Article 

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

Google Scholar  More

Research model and softwareCLIMEX modelFAW growth and development are primarily related to climate conditions, especially temperature patterns17. The current study used CLIMEX (version 4)42, a semi-mechanistic niche modeling platform, to project FAW distribution in relation to climate. The model parameters that describe the species’ response to climate were overlaid onto FAW occurrence data and climate data to project the species’ potential global distribution. Briefly, the annual growth index (GI) was used to describe the potential for FAW population growth during favorable climatic conditions, while stress indices (SI: cold, wet, hot, and dry) and interaction stresses (SX: hot-dry, hot-wet, cold-dry, and cold-wet) (Table 1) were applied to describe the probability that FAW populations could survive unfavorable conditions. The Ecoclimatic index (EI) was derived from a combination of GI, SI, and SX indices to provide an overall annual index of climatic suitability on a scale of 0–10042. An EI value of 0 indicates that the location is not suitable for the long-term survival of the species, whereas an EI value of 100 indicates maximum climatic suitability comparable to conditions in incubators. EI values of more than 30 indicate the optimal climate for a species. In this study, the climatic suitability was classified into four arbitrary categories; unsuitable for EI = 0, marginal for 0  More

1.Stein, E. D., Cohen, Y. & Winer, A. M. Environmental distribution and transformation of mercury compounds. Crit. Rev. Environ. Sci. Technol. 26, 1–43 (1996).CAS 
Article 

Google Scholar 
2.Ciccarelli, C. et al. Assessment of sampling methods about level of mercury in fish. Ital. J. Food Saf. 8, 153–157 (2019).
Google Scholar 
3.Ditri, F. M. Mercury contamination: What we have learned since Minamata. Environ. Monit. Assess. 19, 165–182 (1991).CAS 
Article 

Google Scholar 
4.Monteiro, L. R. & Furness, R. W. Seabirds as monitors of mercury in the marine environment. Water Air Soil Pollut. 80, 851–870 (1995).CAS 
Article 
ADS 

Google Scholar 
5.Pitter, P. In Hydrochemie 5th edn (ed. Pitter, P.) (VSCHT Praha, 2015).
Google Scholar 
6.Hylander, L. D. & Meili, M. 500 years of mercury production: Global annual inventory by region until 2000 and associated emissions. Sci. Total. Environ. 304, 13–27 (2003).CAS 
Article 
ADS 

Google Scholar 
7.Pacyna, E. G. et al. Global emission of mercury to the atmosphere from anthropogenic sources in 2005 and projections to 2020. Atmos. Environ. 44, 2487–2499 (2010).CAS 
Article 
ADS 

Google Scholar 
8.Pai, P., Niemi, D. & Powers, B. A North American inventory of anthropogenic mercury emissions. Fuel Process. Technol. 65, 101–115 (2000).Article 

Google Scholar 
9.Wang, Q. R., Kim, D., Dionysiou, D. D., Sorial, G. A. & Timberlake, D. Sources and remediation for mercury contamination in aquatic systems: A literature review. Environ. Pollut. 131, 323–336 (2004).Article 

Google Scholar 
10.Buck, D. G. et al. A global-scale assessment of fish mercury concentrations and the identification of biological hotspots. Sci. Total Environ. 687, 956–966 (2019).CAS 
Article 
ADS 

Google Scholar 
11.Gentes, S. et al. Application of European water framework directive: Identification reference sites and bioindicator fish species for mercury in tropical freshwater ecosystems (French Guiana). Ecol. Indic. 106, 105468. https://doi.org/10.1016/j.ecolind.2019.105468 (2019).CAS 
Article 

Google Scholar 
12.Thomas, S. M. et al. Climate and landscape conditions indirectly affect fish mercury levels by altering lake water chemistry and fish size. Environ. Res. 188, 109750. https://doi.org/10.1016/j.envres.2020.109750 (2020).CAS 
Article 
PubMed 

Google Scholar 
13.Zupo, V. et al. Mercury accumulation in freshwater and marine fish from the wild and from aquaculture ponds. Environ. Pollut. 255, 112975. https://doi.org/10.1016/j.envpol.2019.112975 (2019).CAS 
Article 
PubMed 

Google Scholar 
14.Zhang, J. L. et al. Health risk assessment of heavy metals in Cyprinus carpio (Cyprinidae) from the upper Mekong river. Environ. Sci. Pollut. Res. 26, 9490–9499 (2019).CAS 
Article 

Google Scholar 
15.Cerna, M. Opatreni mezinarodnich instituci a Ceske republiky k omezovani rizika znecistovani zivotniho prostredi rtuti. Chem. Listy. 98, 916–921 (2004) ((Article in Czech)).CAS 

Google Scholar 
16.Janouskova, D. & Svehla, J. Mercury concentrations in fish tissues in the water reservoir Rimov, South Bohemia. Crop Sci. 19, 43–48 (2002).
Google Scholar 
17.Purba, J. S., Silalahi, J. & Haro, G. Analysis of mercury in fish, North Sumatera, Indonesia by atomic absorption spectrophotometer. Asian J. Pharm. 8, 21–25 (2020).CAS 
Article 

Google Scholar 
18.Willacker, J. J., Eagles-Smith, C. A. & Blazer, V. S. Mercury bioaccumulation in freshwater fishes of the Chesapeake Bay watershed. Ecotoxicology 29, 459484 (2020).Article 

Google Scholar 
19.Regulation (EU) 2017/852 of European Parliament and of the council of 17 May 2017 on mercury, and repealing Regulation (EC) No 1102/2008. Official Journal of the European Union.20.European Commission. The EU Fish Market. https://www.eumofa.eu/documents/20178/415635/EN_The+EU+fish+market_2020.pdf (2020).21.Nebesky, V., Policar, T., Blecha, M., Kristan, J. & Svacina, P. Trends in import and export of fishery products in the Czech Republic during 2010–2015. Aquacult. Int. 24, 1657–1668 (2016).Article 

Google Scholar 
22.FAO. Fisheries & Aquaculture—National Aquaculture Sector Overview—Czech Republic. http://www.fao.org/fishery/countrysector/naso_czechrepublic/en (accessed April 24 April 2021) (2021).23.Rakmanikhah, Z., Esmaili-Sari, A. & Bahramifar, N. Total mercury and methylmercury concentrations in native and invasive fish species in Shadegan International Wetland, Iran, and health risk assessment. Environ. Sci. Pollut. Res. 27, 6765–6773 (2020).Article 

Google Scholar 
24.Celechovska, O., Svobodova, Z., Zlabek, V. & Macharackova, B. Distribution of metals in tissues of the common carp (Cyprinus carpio L.). Acta Vet. Brno 76, 93–100 (2007).Article 

Google Scholar 
25.Cerveny, D. et al. Fish fin-clips as non-lethal approach for biomonitoring of mercury contamination in aquatic environments and human health risk assessment. Chemosphere 163, 290–295 (2016).CAS 
Article 
ADS 

Google Scholar 
26.WHO. Evaluations of the Joint FAO/WHO Expert Committee on Food Additives (JECFA). https://apps.who.int/food-additives-contaminants-jecfa-database/search.aspx.27.Kannan, K. et al. Distribution of total mercury and methyl mercury in water, sediment, and fish from south Florida estuaries. Arch. Environ. Con. Tox. 34, 109–118 (1998).CAS 
Article 

Google Scholar 
28.US EPA. Guidance for Assessing Chemical Contaminant Data for Use in Fish Advisories Documents. Volume 2: Risk Assessment and Fish Consumption Limits, Third Edition. https://www.epa.gov/fish-tech/guidance-assessing-chemical-contaminant-data-use-fish-advisories-documents (accessed 8 May 2021) (2000).29.Ministry of Agriculture of the Czech Republic. Situacni a vyhledova zprava—Ryby. http://eagri.cz/public/web/file/666957/Ryby_2020_web.pdf (accessed 8 May 2021, in Czech) (2020).30.Novotna, K., Svobodova, Z., Harustiakova, D. & Mikula, P. Spatial and temporal trends in contamination of the Czech part of the Elbe River by mercury between 1991 and 2016. Bull. Environ. Contam. Toxicol. 105, 750–757 (2020).CAS 
Article 

Google Scholar 
31.Raldua, D., Diez, S., Bayona, J. M. & Barcelo, D. Mercury levels and liver pathology in feral fish living in the vicinity of a mercury cell chlor-alkali factory. Chemosphere 66, 1217–1225 (2007).CAS 
Article 
ADS 

Google Scholar 
32.Squadrone, S. et al. Heavy metals distribution in muscle, liver, kidney and gill of European catfish (Silurus glanis) from Italian rivers. Chemosphere 90, 358–365 (2013).CAS 
Article 
ADS 

Google Scholar 
33.Cerveny, D. et al. Contamination of fish in important fishing grounds of the Czech Republic. Ecotoxicol. Environ. Saf. 109, 101–109 (2014).CAS 
Article 

Google Scholar 
34.Marsalek, P., Svobodova, Z. & Randak, T. The content of total mercury and methylmercury in common carp from selected Czech ponds. Aquac. Int. 15, 299–304 (2007).CAS 
Article 

Google Scholar 
35.Vicarova, P., Docekalova, H., Ridoskova, A. & Pelcova, P. Heavy metals in the common carp (Cyprinus carpio L.) from three reservoirs in the Czech Republic. Czech J. Food Sci. 34, 422–428 (2016).CAS 
Article 

Google Scholar 
36.Akerblom, S., Bignert, A., Meili, M., Sonesten, L. & Sundbom, M. Half a century of changing mercury levels in Swedish freshwater fish. Ambio 43, 91–103 (2014).Article 

Google Scholar 
37.Dvorak, P., Andreji, J., Mraz, J. & Dvorakova Liskova, Z. Concentration of heavy and toxic metals in fish and sediments from the Morava river basin, Czech Republic. Neuroendocrinol. Lett. 36, 126–132 (2015).CAS 
PubMed 

Google Scholar 
38.Dusek, L. et al. Bioaccumulation of mercury in muscle tissue of fish in the Elbe River (Czech Republic): Maultispecies monitoring study 1991–1996. Ecotoxicol. Environ. Saf. 61, 256–267 (2005).CAS 
Article 

Google Scholar 
39.Marsalek, P., Svobodova, Z. & Randak, T. Total mercury and methylmercury contamination in fish from various sites along the Elbe River. Acta Vet. Brno. 75, 579–585 (2006).CAS 
Article 

Google Scholar 
40.Wang, X. & Wang, W. X. The three ‘B’ of mercury in China: Bioaccumulation, biodynamics and biotransformation. Environ. Pollut. 250, 216–232 (2019).CAS 
Article 

Google Scholar 
41.Jankovska, I. et al. Importance of fish gender as a factor in environmental monitoring of mercury. Environ. Sci. Pollut. Res. 21, 6239–6242 (2014).CAS 
Article 

Google Scholar 
42.Carrasco, L. et al. Patterns of mercury and methylmercury bioaccumulation in fish species downstream of a long-term mercury-contaminated site in the lower Ebro River (NE Spain). Chemosphere 84, 1642–1649 (2011).CAS 
Article 
ADS 

Google Scholar 
43.Havelkova, M., Dusek, L., Nemethova, D., Poleszczuk, G. & Svobodova, Z. Comparison of mercury distribution between liver and muscle: A biomonitoring of fish from lightly and heavily contaminated localities. Sensors. 8, 4095–4109 (2008).CAS 
Article 
ADS 

Google Scholar 
44.Kruzikova, K. et al. The correlation between fish mercury liver/muscle ratio and high and low levels of mercury contamination in Czech localities. Int. J. Electrochem. Sc. 8, 45–56 (2013).CAS 

Google Scholar 
45.Kensova, R., Kruzikova, K. & Svobodova, Z. Mercury speciation and safety of fish from important fishing locations in the Czech Republic. Czech J. Food Sci. 30, 276–284 (2012).CAS 
Article 

Google Scholar 
46.European Commission. Commission Regulation 1881/2006 Setting Maximum Levels of Certain Contaminants in Foodstuffs. https://eur-lex.europa.eu/ (accessed 2 May 2021) (2006). More

Deer capture and samplingWe captured 100 mule deer (54 females, 46 males) during November 2018–February 2019, avoiding capture and sampling of juveniles. We attempted to distribute captures throughout the ~23 km2 study area described by Miller et al.5 to minimize spatial disparities in comparing contemporary and past data, and to assure marks were widely distributed for December ground counts to estimate deer abundance5,35,36. Field and sampling methods generally followed those used elsewhere5,31,37. Field procedures were reviewed and approved by the CPW Animal Care and Use Committee (file 14–2018).We pursued deer on foot and darted them opportunistically, delivering sedative combinations intramuscularly via projectile syringe. Premixed immobilization drug combinations included either nalbuphine (N; 0.9 mg/kg) or butorphanol (B; 0.5 mg/kg) combined with azaperone (A; 0.2 mg/kg) and medetomidine (M; 0.2 mg/kg)38, with standard total doses for respective combinations based on an estimated mass of 70 kg (average drug volume per animal was 1.3 ml NMA, 1.4 ml BAM). We collected rectal mucosa biopsies to determine CWD infection status37. We also collected whole blood and marked all deer with individually identifiable ear tags and some with telemetry (n = 51) or visual identification (n = 12) collars. Ages were estimated to the nearest year via tooth replacement and wear patterns39; observers used a pocket reference guide in the field to help assure consistency. To antagonize sedation upon completion of handling and sampling, each deer received 5 mg atipamezole/mg M administered, injected intramuscularly.Prion diagnosticsFormalin-fixed tissue biopsies were processed and analyzed by immunohistochemistry (IHC) at the Colorado State University Veterinary Diagnostic Laboratory (Fort Collins, Colorado USA; CSUVDL) for evidence of CWD-associated prion (PrPCWD) accumulations using monoclonal antibody F99/97.6.1 (VMRD Inc., Pullman, Washington, USA)40 and standard IHC methods24,37,41, except that the CSUVDL’s IHC staining machine (Leica Microsystems Inc., Buffalo Grove, Illinois, USA) was different from that used in earlier studies (Ventana Medical Systems, Oro Valley, Arizona, USA). Biopsies were evaluated microscopically and classified as positive (infected) or not detected (negative) based on PrPCWD presence or absence; the same pathologist (T. R. Spraker) read biopsies for both the current and prior5 studies.We included only data from deer with biopsies providing ≥3 lymphoid follicles in analyses involving infection status in order to maintain a relatively high (≥90%) probability of detecting infected individuals24. Two animals with low follicle counts that died shortly after capture were excepted by substituting postmortem IHC results. Limiting the acceptable follicle count excluded seven females (two 225SS, five 225SF) and two males (one 225SS, one 225SF) from some analyses. One male deer was 225FF and one female deer was missing a blood sample and thus not assigned to a PRNP gene group; these two individuals also were excluded from some analyses (e.g., Table 1).
PRNP genotypingWe used DNA extracted from whole blood buffy coat aliquots (n = 99) to screen for the presence of sequences at PRNP gene codon 225 that encode for serine (S) and/or phenylalanine (F) in the mature prion polypeptide, classifying individuals as 225SS, 225SF, or 225FF16,36,42. Methods generally followed those described by Jewell et al.16. Briefly, we extracted DNA using the DNeasy® blood and tissue kit (Qiagen, Valenica, California, USA). We amplified the complete open reading frame (ORF) plus 25 bp of 5′ flanking sequences and 53 bp of 3′ flanking sequences in the PRNP coding region using polymerase chain reaction (PCR). Purified DNA was combined in a 0.2 ml PCR tube containing a puReTag Ready-To-Go PCR bead (illustra™, GE Healthcare Bio-Sciences Corp, Piscataway, New Jersey, USA). Each PCR bead contained 2.5 units puReTag DNA polymerase, 10 mM Tris-HCI, 50 mM KCl, 1.5 mM MgCl2, 200 µM of each deoxynucleoside triphosphate, and stabilizers, including bovine serum albumin. For each PCR assay, 1 μL of each primer at 200 nM, 22 μL of RNase-free water and 1 μL of approximately 100 ng total genomic DNA was added for a final volume of 25 μL. Primers used for amplification were forward (MD582F, 5′-ACATGGGCATATGATGCTGACACC-3′) and reverse (MD1479RC, 5′-ACTACAGGGCTGCAGGTAGATACT-3′) described by Jewell et al.16. Reactions were thermal-cycled in a PTC 100 (MJ Research) at 94 C for 5 min and then 32 cycles of 94 C for 7 s, 62 C for 15 s, 72 C for 30 s and a final cycle of 72 C for 5 min, and kept at 4 C until inspected for successful amplification by agarose gel electrophoresis. As confirmed by LaCava et al.19, the MD582F and MD1479RC primers developed by Jewell et al.16 specifically amplify the functional PRNP gene ORF, thereby excluding confounding effects that could arise from the presence of a processed pseudogene that occurs in a majority of deer (Odocoileus spp.)42.We used EcoRI restriction digestion of the PCR-amplified PRNP region16—a validated assay targeting the singular polymorphism at codon 225 in mule deer—to screen all 99 samples for presence of S or F codons. Aliquots (10 μl) of completed PCR reactions were incubated with 10 U EcoRI (New England Biolabs) in a total volume of 12 μl containing 50 mM NaCl, 100 mM Tris/HCl, 10 mM MgCl2, 0.025% Triton X-100 (pH 7.5) at 37 C for 2–16 h followed by the addition of 2.5 μl 6× concentrate gel loading solution (Sigma- Aldrich) per sample, and the inspection of products by agarose gel electrophoresis for the presence of one 897bp-sized band for 225SS, two bands—one 897 bp and one 719 bp—for 225SF, or one 719 bp-sized band for 225FF. As noted by Jewell et al.16, occurrence of TTC (the F codon) at position 225 creates an EcoRI recognition DNA sequence and cleavage site GAATTC from codons 224–225, whereas TCC (the S codon) creates the sequence GAATCC, which is not cut by EcoRI. When incubated with EcoRI, PCR products with a TTC codon at position 225 yielded cleavage fragments of the predictable sizes listed16. Because no other sites within the PRNP ORF DNA sequence are potentially transformable to GAATTC with one base change, this represents a specific genotyping method for assessing the S225F polymorphism in mule deer16.To confirm findings from EcoRI screening, we examined sequences of the complete PRNP ORF from 20 samples that showed evidence of cleavage indicating 225*F and 6 samples without cleavage identified as 225SS. For DNA sequencing, we used primers 245 (5′-GGTGGTGACTGACTGTGTGTTGCTTGA-3′), 12 (5′-TGGTGGTGACTGTGTGTTGCTTGA-3′) and 3FL1 (5′-GATTAAGAAGATAATGAAAACAGGAAGG-3′; Integrated DNA Technologies). Sanger sequencing was done on purified PCR product by Eurofins Genomics (Louisville, Kentucky, USA). Sequence chromatograms were viewed and DNA sequence alignments and comparisons were made using the MAFFT multiple sequence alignment program v7.450 module, software platform v2020.2.3 of Geneious Prime. Sequencing confirmed the presence of coding for F in all samples identified as 225*F by EcoRI digestion, as well as the absence of such coding in samples identified by EcoRI digestion as 225SS. Moreover, presence of AGC at codon 138 in all sequenced samples reconfirmed that the primers we used had amplified the functional PRNP gene42.Statistics and reproducibilityFor analyses, we tabulated IHC-positive and -negative results to estimate apparent prevalence of prion infection. We also tabulated the number of individuals assigned to PRNP genotypes and to age groupings as described. Age groupings were selected based on relevance to CWD epidemiology in mule deer1,5,8,12,16,17,18,20,24,31,37. Assuming a ~2-year disease course5,8,17 and relative scarcity of end-stage disease in 225SS deer More