Ecology

D’Costa, V. M. et al. Antibiotic resistance is ancient. Nature 477, 457–461. https://doi.org/10.1038/nature10388 (2011).Article 
ADS 

Google Scholar 
Cytryn, E. The soil resistome: The anthropogenic, the native, and the unknown. Soil Biol. Biochem. 63, 18–23. https://doi.org/10.1016/j.soilbio.2013.03.017 (2013).Article 

Google Scholar 
Holmes, A. H. et al. Understanding the mechanisms and drivers of antimicrobial resistance. Lancet 387, 176–187. https://doi.org/10.1016/S0140-6736(15)00473-0 (2016).Article 

Google Scholar 
Regulation (EC) No 1831/2003 of the European parliament and of the council of 22 September 2003 on additives for use in animal nutrition.European Commission. Communication from the commission to the European parliament, the council, the European economic and social committee and the committee of the regions: A farm to fork strategy for a fair, healthy and environmentally-friendly food system COM/2020/381 Final, (2020).Kumar, K. C., Gupta, S. C., Chander, Y. & Singh, A. K. Antibiotic use in agriculture and its impact on the terrestrial environment. Adv. Agron. 87, 1–54. https://doi.org/10.1016/S0065-2113(05)87001-4 (2005).Article 

Google Scholar 
Chee-Sanford, J. C. et al. Fate and transport of antibiotic residues and antibiotic resistance genes following land application of manure waste. J. Environ. Qual. 38, 1086–1108. https://doi.org/10.2134/jeq2008.0128 (2009).Article 

Google Scholar 
Heuer, H., Schmitt, H. & Smalla, K. Antibiotic resistance gene spread due to manure application on agricultural fields. Curr. Opin. Microbiol. 14, 236–243. https://doi.org/10.1016/j.mib.2011.04.009 (2011).Article 

Google Scholar 
Epelde, L. et al. Characterization of composted organic amendments for agricultural use. Front. Sustain. Food Syst. 2, 44. https://doi.org/10.3389/fsufs.2018.00044 (2018).Article 

Google Scholar 
Youngquist, C. P., Mitchell, S. M. & Cogger, C. G. Fate of antibiotics and antibiotic resistance during digestion and composting: A review. J. Environ. Qual. 45, 537–545. https://doi.org/10.2134/jeq2015.05.0256 (2016).Article 

Google Scholar 
Ma, X., Xue, X., González-Mejía, A., Garland, J. & Cashdollar, J. Sustainable water systems for the city of tomorrow: A conceptual framework. Sustainability 7, 12071–12105. https://doi.org/10.3390/su70912071 (2015).Article 

Google Scholar 
Wang, Y. et al. Degradation of antibiotic resistance genes and mobile gene elements in dairy manure anerobic digestion. PLoS ONE 16, e0254836. https://doi.org/10.1371/journal.pone.0254836 (2021).Article 

Google Scholar 
Thanomsub, B. et al. Effects of ozone treatment on cell growth and ultrastructural changes in bacteria. J. Gen. Appl. Microbiol. 48, 193–199. https://doi.org/10.2323/jgam.48.193 (2002).Article 

Google Scholar 
Sousa, J. M. et al. Ozonation and UV254nm radiation for the removal of microorganisms and antibiotic resistance genes from urban wastewater. J. Hazard. Mater. 323, 434–441. https://doi.org/10.1016/j.jhazmat.2016.03.096 (2017).Article 

Google Scholar 
Park, J. H., Choppala, G. K., Bolan, N. S., Chung, J. W. & Chuasavathi, T. Biochar reduces the bioavailability and phytotoxicity of heavy metals. Plant Soil 348, 439–451. https://doi.org/10.1007/s11104-011-0948-y (2011).Article 

Google Scholar 
Jeffery, S. et al. The way forward in biochar research: targeting trade-offs between the potential wins. GCB Bioenergy 7, 1–13. https://doi.org/10.1111/gcbb.12132 (2015).Article 

Google Scholar 
Krasucka, P. et al. Engineered biochar: A sustainable solution for the removal of antibiotics from water. Chem. Eng. J. 405, 126926. https://doi.org/10.1016/j.cej.2020.126926 (2021).Article 

Google Scholar 
Ken, D. S. & Sinha, A. Recent developments in surface modification of Nano zero-valent iron (nZVI): remediation, toxicity and environmental impacts. Environ. Nanotechnol. Monit. Manag. 14, 100344. https://doi.org/10.1016/j.enmm.2020.100344 (2020).Article 

Google Scholar 
Zhao, X. et al. An overview of preparation and applications of stabilized zero-valent iron nanoparticles for soil and groundwater remediation. Water Res. 100, 245–266. https://doi.org/10.1016/j.watres.2016.05.019 (2016).Article 

Google Scholar 
Diao, M. & Yao, M. Use of zero-valent iron nanoparticles in inactivating microbes. Water Res. 43, 5243–5251. https://doi.org/10.1016/j.watres.2009.08.051 (2009).Article 

Google Scholar 
Shi, C. J., Wei, J., Jin, Y., Kniel, K. E. & Chiu, P. C. Removal of viruses and bacteriophages from drinking water using zero-valent iron. Sep. Purif. Technol. 84, 72–78. https://doi.org/10.1016/j.seppur.2011.06.036 (2012).Article 

Google Scholar 
Anza, M., Salazar, O., Epelde, L., Alkorta, I. & Garbisu, C. The application of nanoscale zero-valent iron promotes soil remediation while negatively affecting soil microbial biomass and activity. Front. Environ. Sci. 7, 19. https://doi.org/10.3389/fenvs.2019.00019 (2019).Article 

Google Scholar 
FAOSTAT. Mushrooms and truffles, production quantity (tons). https://www.tilasto.com/en/topic/geography-and-agriculture/crop/mushrooms-and-truffles/mushrooms-and-truffles-production-quantity/spain, (2020).Polat, E., Uzun, H., Topçuo, B., Önal, K. & Onus, A. N. Effects of spent mushroom compost on quality and productivity of cucumber (Cucumis sativus L.) grown in greenhouses. Afr. J. Biotechnol. 8, 176–180 (2009).
Google Scholar 
Fazaeli, H. & Masoodi, A. R. T. Spent wheat straw compost of Agaricus bisporus mushroom as ruminant feed. Asian-Australas. J. Anim. Sci. 19, 845–851. https://doi.org/10.5713/ajas.2006.845 (2006).Article 

Google Scholar 
Phan, C. W. & Sabaratnam, V. Potential uses of spent mushroom substrate and its associated lignocellulosic enzymes. Appl. Microbiol. Biotechnol. 96, 863–873. https://doi.org/10.1007/s00253-012-4446-9 (2012).Article 

Google Scholar 
Lau, K. L., Tsang, Y. Y. & Chiu, S. W. Use of spent mushroom compost to bioremediate PAH-contaminated samples. Chemosphere 52, 1539–1546. https://doi.org/10.1016/S0045-6535(03)00493-4 (2003).Article 
ADS 

Google Scholar 
Mayans, B. et al. An assessment of Pleurotus ostreatus to remove sulfonamides, and its role as a biofilter based on its own spent mushroom substrate. Environ. Sci. Pollut. Res. Int. 28, 7032–7042. https://doi.org/10.1007/s11356-020-11078-3 (2021).Article 

Google Scholar 
Congilosi, J. L. & Aga, D. S. Review on the fate of antimicrobials, antimicrobial resistance genes, and other micropollutants in manure during enhanced anaerobic digestion and composting. J. Hazard. Mater. 405, 123634. https://doi.org/10.1016/j.jhazmat.2020.123634 (2021).Article 

Google Scholar 
Oliver, J. P. et al. Invited review: fate of antibiotic residues, antibiotic-resistant bacteria, and antibiotic resistance genes in US dairy manure management systems. J. Dairy Sci. 103, 1051–1071. https://doi.org/10.3168/jds.2019-16778 (2020).Article 

Google Scholar 
Beneragama, N. et al. Survival of multidrug-resistant bacteria in thermophilic and mesophilic anaerobic co-digestion of dairy manure and waste milk. Anim. Sci. J. 84, 426–433. https://doi.org/10.1111/asj.12017 (2013).Article 

Google Scholar 
Sun, W., Qian, X., Gu, J., Wang, X. J. & Duan, M. L. Mechanism and effect of temperature on variations in antibiotic resistance genes during anaerobic digestion of dairy manure. Sci. Rep. 6, 30237. https://doi.org/10.1038/srep30237 (2016).Article 
ADS 

Google Scholar 
Sun, W., Gu, J., Wang, X., Qian, X. & Peng, H. Solid-state anaerobic digestion facilitates the removal of antibiotic resistance genes and mobile genetic elements from cattle manure. Bioresour. Technol. 274, 287–295. https://doi.org/10.1016/j.biortech.2018.09.013 (2019).Article 

Google Scholar 
Zou, Y., Xiao, Y., Wang, H., Fang, T. & Dong, P. New insight into fates of sulfonamide and tetracycline resistance genes and resistant bacteria during anaerobic digestion of manure at thermophilic and mesophilic temperatures. J. Hazard. Mater. 384, 121433. https://doi.org/10.1016/j.jhazmat.2019.121433 (2020).Article 

Google Scholar 
Agga, G. E., Kasumba, J., Loughrin, J. H. & Conte, E. D. Anaerobic digestion of tetracycline spiked livestock manure and poultry litter increased the abundances of antibiotic and heavy metal resistance genes. Front Microbiol. 11, 614424. https://doi.org/10.3389/fmicb.2020.614424 (2020).Article 

Google Scholar 
Jauregi, L., Epelde, L., González, A., Lavín, J. L. & Garbisu, C. Reduction of the resistome risk from cow slurry and manure microbiomes to soil and vegetable microbiomes. Environ. Microbiol. 23, 7643–7660. https://doi.org/10.1111/1462-2920.15842 (2021).Article 

Google Scholar 
Zhang, Z. et al. Assessment of global health risk of antibiotic resistance genes. Nat Commun 13, 1553. https://doi.org/10.1038/s41467-022-29283-8 (2022).Article 
ADS 

Google Scholar 
He, Y. et al. Antibiotic resistance genes from livestock waste: occurrence, dissemination, and treatment. npj Clean Water 3, 4. https://doi.org/10.1038/s41545-020-0051-0 (2020).Article 

Google Scholar 
Cui, E., Wu, Y., Zuo, Y. & Chen, H. Effect of different biochars on antibiotic resistance genes and bacterial community during chicken manure composting. Bioresour. Technol. 203, 11–17. https://doi.org/10.1016/j.biortech.2015.12.030 (2016).Article 

Google Scholar 
Fu, Y., Zhang, A., Guo, T., Zhu, Y. & Shao, Y. Biochar and hyperthermophiles as additives accelerate the removal of antibiotic resistance genes and mobile genetic elements during composting. Materials (Basel) 14, 5428. https://doi.org/10.3390/ma14185428 (2021).Article 
ADS 

Google Scholar 
Forsberg, K. J. et al. Bacterial phylogeny structures soil resistomes across habitats. Nature 509, 612–616. https://doi.org/10.1038/nature13377 (2014).Article 
ADS 

Google Scholar 
Li, H. et al. Effects of bamboo charcoal on antibiotic resistance genes during chicken manure composting. Ecotoxicol. Environ. Saf. 140, 1–6. https://doi.org/10.1016/j.ecoenv.2017.01.007 (2017).Article 
ADS 

Google Scholar 
Bondarenko, O., Ivask, A., Käkinen, A. & Kahru, A. Sub-toxic effects of CuO nanoparticles on bacteria: Kinetics, role of Cu ions and possible mechanisms of action. Environ. Pollut. 169, 81–89. https://doi.org/10.1016/j.envpol.2012.05.009 (2012).Article 

Google Scholar 
Wang, X. et al. Bacterial exposure to ZnO nanoparticles facilitates horizontal transfer of antibiotic resistance genes. NanoImpact 10, 61–67. https://doi.org/10.1016/j.impact.2017.11.006 (2018).Article 
ADS 

Google Scholar 
Qiu, X., Zhou, G. & Wang, H. Nanoscale zero-valent iron inhibits the horizontal gene transfer of antibiotic resistance genes in chicken manure compost. J. Hazard. Mater. 422, 126883. https://doi.org/10.1016/j.jhazmat.2021.126883 (2022).Article 

Google Scholar 
Zeng, T., Wilson, C. J. & Mitch, W. A. Effect of chemical oxidation on the sorption tendency of dissolved organic matter to a model hydrophobic surface. Environ. Sci. Technol. 48, 5118–5126. https://doi.org/10.1021/es405257b (2014).Article 
ADS 

Google Scholar 
Pak, G. et al. Comparison of antibiotic resistance removal efficiencies using ozone disinfection under different pH and suspended solids and humic substance concentrations. Environ. Sci. Technol. 50, 7590–7600. https://doi.org/10.1021/acs.est.6b01340 (2016).Article 
ADS 

Google Scholar 
Zhuang, Y. et al. Inactivation of antibiotic resistance genes in municipal wastewater by chlorination, ultraviolet, and ozonation disinfection. Environ. Sci. Pollut. Res. Int. 22, 7037–7044. https://doi.org/10.1007/s11356-014-3919-z (2015).Article 

Google Scholar 
Park, S., Rana, A., Sung, W. & Munir, M. Competitiveness of quantitative polymerase chain reaction (qPCR) and droplet digital polymerase chain reaction (ddPCR) technologies, with a particular focus on detection of antibiotic resistance genes (ARGs). Appl. Microbiol. 1, 426–444. https://doi.org/10.3390/applmicrobiol1030028 (2021).Article 

Google Scholar 
European Medicines Agency. European surveillance of veterinary antimicrobial consumption, (2020). Sales of Veterinary Antimicrobial Agents in 31 European Countries in 2018 (EMA/24309/2020).Heuer, H. et al. The complete sequences of plasmids pB2 and pB3 provide evidence for a recent ancestor of the IncP-1β group without any accessory genes. Microbiology (Reading) 150, 3591–3599. https://doi.org/10.1099/mic.0.27304-0 (2004).Article 

Google Scholar 
World Health Organization. Critically Important Antimicrobials for Human Medicine, 6th Revision (WHO, Geneva, Switzerland, 2019).Zhu, Y. G. et al. Diverse and abundant antibiotic resistance genes in Chinese swine farms. Proc. Natl Acad. Sci. U. S. A. 110, 3435–3440. https://doi.org/10.1073/pnas.1222743110 (2013).Article 
ADS 

Google Scholar 
Guo, T. et al. Increased occurrence of heavy metals, antibiotics and resistance genes in surface soil after long-term application of manure. Sci. Total Environ. 635, 995–1003. https://doi.org/10.1016/j.scitotenv.2018.04.194 (2018).Article 
ADS 

Google Scholar 
Nõlvak, H. et al. Inorganic and organic fertilizers impact the abundance and proportion of antibiotic resistance and integron-integrase genes in agricultural grassland soil. Sci. Total Environ. 562, 678–689. https://doi.org/10.1016/j.scitotenv.2016.04.035 (2016).Article 
ADS 

Google Scholar 
Chen, Q. L. et al. Effect of biochar amendment on the alleviation of antibiotic resistance in soil and phyllosphere of Brassica chinensis L.. Soil Biol. Biochem. 119, 74–82. https://doi.org/10.1016/j.soilbio.2018.01.015 (2018).Article 

Google Scholar 
Zhu, B., Chen, Q., Chen, S. & Zhu, Y. G. Does organically produced lettuce harbor higher abundance of antibiotic resistance genes than conventionally produced?. Environ. Int. 98, 152–159. https://doi.org/10.1016/j.envint.2016.11.001 (2017) .Article 

Google Scholar 
Métodos, M. A. P. A. Oficiales de análisis de suelos y Aguas Para riego. Minist. Agric. Pesca Aliment. Métodos Oficiales Anal. III (1994).Muziasari, W. I. et al. Aquaculture changes the profile of antibiotic resistance and mobile genetic element associated genes in Baltic Sea sediments. FEMS Microbiol. Ecol. 92, fiw052. https://doi.org/10.1093/femsec/fiw052 (2016).Article 

Google Scholar 
Muurinen, J. et al. Influence of manure application on the environmental resistome under Finnish agricultural practice with restricted antibiotic use. Environ. Sci. Technol. 51, 5989–5999. https://doi.org/10.1021/acs.est.7b00551 (2017).Article 
ADS 

Google Scholar 
Muziasari, W. I. et al. The resistome of farmed fish feces contributes to the enrichment of antibiotic resistance genes in sediments below Baltic Sea fish farms. Front. Microbiol. 7, 2137. https://doi.org/10.3389/fmicb.2016.02137 (2017).Article 

Google Scholar 
Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods 25, 402–408. https://doi.org/10.1006/meth.2001.1262 (2001).Article 

Google Scholar 
Ovreås, L., Forney, L., Daae, F. L. & Torsvik, V. Distribution of bacterioplankton in meromictic Lake Saelenvannet, as determined by denaturing gradient gel electrophoresis of PCR-amplified gene fragments coding for 16S rRNA. Appl. Environ. Microbiol. 63, 3367–3373. https://doi.org/10.1128/aem.63.9.3367-3373.1997 (1997) .Article 
ADS 

Google Scholar 
Caporaso, J. G. et al. Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms. ISME J. 6, 1621–1624. https://doi.org/10.1038/ismej.2012.8 (2012).Article 

Google Scholar 
Lanzén, A. et al. Multi-targeted metagenetic analysis of the influence of climate and environmental parameters on soil microbial communities along an elevational gradient. Sci. Rep. 6, 28257. https://doi.org/10.1038/srep28257 (2016).Article 
ADS 

Google Scholar 
Pinna, N. K., Dutta, A., Monzoorul, H. M. & Mande, S. S. Can targeting non-contiguous V-regions with paired-end sequencing improve 16S rRNA-based taxonomic resolution of microbiomes?: An in silico evaluation. Front. Genet. 10, 653. https://doi.org/10.3389/fgene.2019.00653 (2019).Article 

Google Scholar 
Andrews, S. FastQC: A quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (2010).Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12. https://doi.org/10.14806/ej.17.1.200 (2011).Article 

Google Scholar 
Bolyen, E. et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857. https://doi.org/10.1038/s41587-019-0209-9 (2019).Article 

Google Scholar 
Amir, A. et al. Deblur rapidly resolves single-nucleotide community sequence patterns. mSystems 2, e00191-e216. https://doi.org/10.1128/mSystems.00191-16 (2017).Article 

Google Scholar 
Yang, Y., Li, B., Zou, S., Fang, H. H. P. & Zhang, T. Fate of antibiotic resistance genes in sewage treatment plant revealed by metagenomic approach. Water Res. 62, 97–106. https://doi.org/10.1016/j.watres.2014.05.019 (2014).Article 

Google Scholar 
Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).Book 
MATH 

Google Scholar 
de Mendiburu, F. Agricolae: Statistical procedures for agricultural research. R package version 1.3-3. https://CRAN.R-project.org/package=agricolae (2020).Paradis, E. & Schliep, K. ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528. https://doi.org/10.1093/bioinformatics/bty633 (2019).Article 

Google Scholar 
Oksanen, J. et al. Vegan: Community ecology package. R Package Version 2.3-1. (2015). More

Classification of snow algae in the ice core based on ITS2 sequencesWe used high-throughput sequencing to obtain DNA sequences of algae from 19 layers of an ice core drilled on a glacier in central Asia, dated from present time to 8000 years ago (Fig. 1 and Table S1). In total, 17,016 unique sequences (phylotypes) for the fast-evolving algal nuclear rDNA internal transcribed spacer 2 (ITS2) region were determined in the ice core, from which 290 OTUs were defined with ≥98% nt sequence identity among all OTUs.The ITS2 sequences were classified at the species level according to the genetic species concept based on secondary structural differences in the ITS2 region, which correlate with the boundaries of most biological species [38]. The ITS2 sequences from ice core samples were classified into 24 subgroups consisting of 17 chlorophycean, 5 trebouxiophycean, and 2 ulvophycean groups based on their secondary structures and BLASTn results (Fig. S1 and Tables S3–S4). The 17 subgroups of Chlorophyceae were subdivided into 10 subgroups of the Chloromonadinia clade, 1 subgroup of the Monadinia clade (recently treated as the genus Microglena [54]), 3 subgroups of the Reinhardtinia clade, 2 subgroups of the Stephanosphaerinia clade, and 1 subgroup corresponding to an unnamed group (which is related to Ploeotila sp. CCCryo 086-99) (for the clade names, see [55]). Although the Chloromonadinia clade contains several snow species belonging to Chloromonas or Chlainomonas, the 10 subgroups of the Chloromonadinia clade were considered to be Chloromonas. The 5 trebouxiophycean subgroups were composed of 2 subgroups of the Chlorella group, 1 subgroup of the Raphidonema group, 1 subgroup of the Trebouxia group, and 1 subgroup of the Neocystis group. The 2 subgroups of Ulvophyceae were closely related to the genus Chamaetrichon and Planophila, respectively. It is noted that Sanguina (‘Chlamydomonas’-snow group B [6]), Ancylonema, and Mesotaenium, which are snow algal genera found throughout the world [56, 57], were not detected in the ice core samples (Tables S3–S4).Global distribution of the Raphidonema groupTo understand the process by which snow algae form geographically specific population structures and how they migrate globally across the glaciers and snow fields, it is necessary to focus on the microbial species that inhabit the global cryosphere. Previous work elucidated that the Raphidonema group and ‘Chlamydomonas’-snow group B (Sanguina) are the cosmopolitans at both poles [6], but the latter was not detected in ice core samples examined in this study. Therefore, to elucidate the evolutionary history of the Raphidonema group, we further analyzed the ITS2 sequences obtained from the ice core sample as well as the glacier-surface samples from both poles [6] and from the mid-latitudes (samples from 10 sites, obtained in this study) (hereafter, surface samples; Table S2). Members of the Raphidonema group were detected in the older (deep core) layers of the ice core and at the glacier surface of central Asia (Fig. S1 and Tables S3–S4), as well as in the red snow samples from both poles [6]. In central Asia, the Raphidonema group was found in the Russian, Chinese, and Kyrgyz samples but was not detected in the Japanese and Tajik samples (Fig. S1 and Tables S3–S4). Combining these sequences yielded 893,649 reads and 22,389 unique sequences for subsequent detailed analysis (Tables S5–S6). The taxonomic composition of the Raphidonema communities differed among the mid-latitude, ice core, Arctic, and Antarctic samples as determined by PERMANOVA (Table S7). Most of the unique sequences in the Raphidonema group were consistent with an endemic distribution (Tables S8–S10). An average of 77% of the unique sequences of the Raphidonema group were endemic to a specific region (mid-latitude, 96%; Antarctic, 66%; Arctic, 79%), accounting for 40% of the total sequencing reads (mid-latitude, 77%; Antarctic, 74%; Arctic, 22%) (Fig. 2a, b and Tables S9–S10). This result suggested that most of the unique sequences are endemic, indicating that their dispersal has been limited to their respective regions [58,59,60,61].Fig. 2: Distribution types of the Raphidonema group obtained from each region and the ice core based on ITS2 unique sequences.Proportions of unique sequence and number of sequencing reads are shown. a Unique sequences from surface snow and ice-core samples. b Number of sequencing reads from surface snow and ice core samples. c Unique sequences from the indicated locations within the ice core. d Number of sequencing reads of the unique sequences from the indicated locations within the ice core.Full size imageNext, we analyzed the global distribution of the cosmopolitan phylotypes of the Raphidonema group, because a previous study analyzed their distribution only at the poles [6]. Only a limited number of unique sequences were distributed in all regions (mid-latitude, 1.4%; Antarctic, 5.6%; Arctic, 3.1%), accounting for a large proportion of the sequencing reads in polar regions but for only a small proportion in the mid-latitudes (mid-latitude, 2.8%; Antarctic, 20%; Arctic, 55%) (Figs. 2a, b, S2–S3, and Tables S9–S10). The distribution types of the Raphidonema group obtained from each region and the ice core were similar between the USEARCH and DADA2 analyses (Figs. 2, S4). In addition, we note that in ancient samples, post-mortem nt substitutions, such as cytosine to thymine, accumulate over many years of deposition [62], and these are not included in the DADA2 error model, which leads to the elimination of minor sequences in the DADA2 analysis. Therefore, we based our analysis on the results of the USEARCH unique sequences. These results suggested that only a few snow algae in the Raphidonema group were detected in samples from the mid-latitude regions.Snow algae of the Raphidonema group were detected in different ice core layers, corresponding to different time periods. The ice core records revealed that the distribution types of the Raphidonema group have not changed significantly for the last 8000 years, with p = 0.1924 based on a PERMANOVA between the newer (1800–2001 AD) and the older (6000–8000 years before present) layers (Fig. 2c, d). In ice core samples, 77% of the unique sequences of the Raphidonema group were detected only in the ice core samples, accounting for 23% of the total sequencing reads (Fig. S5). Although some of these unique sequences may be artifacts of the post-mortem nt substitution or sequencing errors, because we conducted sequence quality filtering and removed the majority of artifact sequences by removing the singleton clusters, most of the unique sequences in the ice core are not likely to be artifacts, but they could represent endemic phylotypes (Figs. 2a, b, S5).The cosmopolitan phylotypes were detected over a broad period as represented by ice core samples. They were present in approximately similar ratios in the newer and older layers (Fig. 2c, d). The cosmopolitan phylotypes were relatively abundant in the ice core samples (average, 4.0%; range, 0.2–13%), accounting for 13% (0.9–81% in the samples) of the total sequencing reads (Figs. 2c, d and S5).Microevolution of cosmopolitan and endemic phylotypesWe analyzed the evolutionary relationship between cosmopolitan and endemic phylotypes of the Raphidonema group among all snow surface and ice core samples. In total, 22,389 unique sequences of the Raphidonema group were clustered into 170 OTUs that were defined with ≥98% nt sequence identity among sequences within OTUs. The OTU sequences were subdivided into five subgroups (Groups A–E) based on phylogenetic analysis (Figs. S6–S11 and Tables S11–S12). Based on a previous study [63], Groups A–C and Group E were assigned to R. sempervirens and R. nivale, respectively, but Group D was not consistent with any species examined in that study (Fig. S6).The phylotypes were categorized into three subsets: the cosmopolitan phylotypes found at both poles and the mid-latitude regions; the multi-region phylotypes found in any two of the Antarctic, Arctic, and mid-latitude regions; and the endemic phylotypes found in only one of the three regions. Cosmopolitan phylotypes were found in Groups A, B, and C and accounted for 64.6% of the unique sequences. We then analyzed the dispersal of the three groups in detail.MJ networks [47] for the ITS2 sequences in each subgroup revealed that the cosmopolitan phylotypes were located at the center of the networks in Groups A and C that contained any types (endemics, multi-regions, and cosmopolitans) of the phylotypes, whereas the endemic phylotypes were considered to be derived from the cosmopolitan phylotypes (Figs. 3 and S12–S13). Moreover, the outgroup phylotypes were directly connected to the cosmopolitan phylotypes. These findings clearly showed that the cosmopolitan phylotypes were ancestral, whereas the endemic phylotypes were derived. In contrast, there were remarkable differences in the shape of the networks between Group B and the others (Groups A and C). In Group B, the Antarctic endemic phylotypes formed a distinct clade, and multi-region phylotypes seemed to be recently derived from this clade. In addition, the Arctic endemic phylotypes formed another distinct clade. These two Group-B clades split directly from a cosmopolitan phylotype (5.3% of the total sequencing reads). For Groups A and C, however, major portions of the total sequencing reads belonged to cosmopolitan phylotypes in Groups A (48.2%) and C (62.4%), and the endemic and multi-region phylotypes were directly connected to these major cosmopolitan phylotypes in a radial manner—the so-called “star-like” pattern [64]. These contrasting network shapes seem to have been formed as a consequence of the unique evolution of each of these groups. We also found that sequences from ice cores did not represent a basal position (Figs. 3 and S12–S13). This is because the haplotypes found in the modern samples have existed from times earlier than the ice core ages, due to the very small mutation numbers expected to have occurred since the ice core ages. Therefore, detected ice core ages were not included in the molecular evolution calculations of our demographic model. However, the phylogenetic networks themselves do not provide information on the evolutionary time scale. Hence, the ice core samples provide further direct evidence that Raphidonema, especially cosmopolitans belonging to this genus, persistently grew on snow and ice at least during the Holocene, and their ITS2 sequences have not changed over the last 8000 years.Fig. 3: Phylogenetic relationships among phylotypes of the Raphidonema groups.Phylotype networks for ITS2 sequences within Groups A (a), B (b), and C (c) of the Raphidonema group that include the cosmopolitan phylotypes in this study. The median-joining method was used. Circles indicate phylotypes; the size of each circle is proportional to the number of unique sequences. Each notch on the edges represents a mutation. Phylotypes are colored according to geographic region. The arrow represents the phylotype in the outgroup (see Fig. S6).Full size imageReferring to “ancestral” phylotypes as those having a longer history than other, more recently derived phylotypes, it is possible that individuals not closely related can share the same ancestral phylotype. In such cases, if genetically far-related individuals from various geographical regions share the same ancestral phylotype, they appear to be “cosmopolitan” (Fig. S14a). In order to distinguish between these “apparent cosmopolitans”, and “true cosmopolitans” that migrate globally, it is necessary to show that the cosmopolitan and endemic phylotypes have distinct demographic histories rather than being part of a continuous population sharing certain demographic dynamics (Fig. S14). Because phylotype networks are not useful for quantifying the rate(s) of microevolution, we used the coalescent model to quantify phylotype demographics [65]. As numerous phylotypes must be analyzed with this approach, we concentrated on statistical inference based on pairwise comparisons of phylotypes, for which the likelihood can be determined in a practical manner (see Materials and Methods). Histograms for the number of mismatched sites between two phylotypes chosen from a set of phylotypes, which will be called the pairwise mismatch distribution, are shown in Figs. 4 and S15. For Groups A and C, the distribution among cosmopolitans, multi-regions, and endemics was unimodal, in which the modes align from left to right with the order cosmopolitans, multi-regions, and endemics. Rogers and Harpending [48] noted that this “wave” propagation results from the expansion in size of a population, which leads to large mismatches, and the mode shifts to the right (see Fig. 2 of [48]). As time passes, the mode shifts to the left and eventually returns to the origin, i.e., representing a population that has not undergone an expansion event. Rogers and Harpending obtained an approximate solution for the wave and fitted the solution to human mitochondrial sequence data. We improved upon their method based on the coalescent model (see Materials and Methods) and applied it to the ITS2 sequence data for snow algae.Fig. 4: Mismatch distribution based on the number of pairwise differences in each distribution type in Raphidonema groups.The lines represent the observed number of pairwise differences in each distribution type (cosmopolitan, multi-region, endemic) within the Raphidonema Groups A (a), B (b) and C (c). Calculations were performed for all distribution types of Raphidonema Groups A and C, for which various cosmopolitan phylotypes were detected. On the other hand, calculations for only multi-region and endemic phylotypes were performed for Raphidonema group B, because no variation was found in cosmopolitan phylotypes.Full size imageFor Group A, when we fit the single demographic model to all phylotypes, the log-likelihood was –414,487. In contrast, when we fit the demographic model to each subset, that is, cosmopolitans, multi-regions, and endemics, separately, the log-likelihood was –341,964. Because the latter is larger than the former, we fit the model to each subset of phylotypes separately. For Group C, when we fit the demographic model to the cosmopolitans, multi-regions, and endemics separately, the log-likelihood was –142,106, which is larger than the log-likelihood, –218,080, when we fit the single demographic model to all phylotypes. In contrast to Groups A and C then, we fit the single demographic model to all phylotypes of Group B because the log-likelihood, –196,070, was larger than the log-likelihood, –220,145, when we fit the demographic model to the cosmopolitans, multi-regions, and endemics separately. These results suggested that cosmopolitans, multi-regions, and endemics experienced different demographic histories in Groups A and C, whereas they had the same demographic history in Group B (Table S13). These results indicate the cosmopolitans in Group A and C are true cosmopolitans, whereas the those in Group B can be regarded as an apparent cosmopolitan.The ML estimates of (tau = 2ut_0), (theta _0 = 2N_0u), and (theta _1 = 2N_1u) are shown in Table S13 with standard deviation values. The population expanded t years ago, with the size before and after the expansion being represented by N0 and N1, respectively. The mutation rate (u) was assumed to be 7.9 × 10–8/ sequence/generation, and the generation interval was assumed to be 24 days (Materials and Methods). In Group A, for the cosmopolitans, the estimates of t, N0, and N1 were (33.8/(2 times 7.9) times 10^8 times {textstyle{{24} over {365}}} = 1.4 times 10^7) years, ((0.108 – 0.010)/(2 times 7.9) times 10^8 = (6.8 – 0.63) times 10^5), and ((0.217)/(2 times 7.9) times 10^8 = 1.4 times 10^6), respectively. In the same way, we computed estimates of t, N0, and N1 of other phylotypes and other groups (Table S14). For the endemics, the respective values were 9.2 × 106 years, 80, and 2.1 × 107, and the values were 4.6 × 106 years, 139, and 1.5 × 107 for the multi-regions. Taking into account the minimum and maximum ranges of the mutation rates per generation as well as the generation intervals, t for cosmopolitans was 3.6 × 106–4.0 × 107 years ago, and t for endemics was 2.3 × 106–2.6 × 107 years ago (Table S14). These results suggested that the cosmopolitans existed at least 1.4 × 107 years ago, and the endemics were derived from the cosmopolitans 9.2 × 106 years ago. The size of the endemics expanded 2.6 × 105-fold, which may have resulted from extensive dispersal. The multi-regions tended to mimic the endemics. Note that our demographic model was simplified to avoid overparameterization. In reality, considering the branching patterns of the MJ network, it is plausible that the endemic phylotypes have been repetitively and continuously derived from the cosmopolitans in multiple lineages—from 9.2 × 106 years ago to the present. In the same way, as for Group C, our results suggested that the cosmopolitan population expanded 3.9-fold ~3.2 × 106 years ago, and the endemics were derived from the cosmopolitans 1.9 × 105 years ago. The size of the endemics expanded 59-fold. In contrast to the phylotypes of Groups A and C, those of Group B experienced no significant expansion (Supplementary Results). In Groups A and C, the derived endemics (and multi-regions) expanded greatly as compared with the ancestral cosmopolitans (Table S14). These extraordinary expansions constitute evidence for local adaptation by the endemic/multi-region populations. In contrast, there was no evidence of local adaptation in Group B. The mismatch distribution of the entire Group B (multi-regions + endemics) showed a multimodal pattern (Fig. 4), which is present in the populations with stable sizes for a long period. When the populations finally reach equilibrium, the mismatch distributions show the exponential distribution [48]. Based on our ML estimates (Table S14), the historical population of Group B has been stable. More

The Campbell et al. model provides an excellent framework to identify species whose extinction leads to community collapse and species whose reintroduction can restore the community (see Fig. 2 for an illustration of these processes). Our first objective, finding the effect of species extinction on the rest of the species in an established community, is achievable using the concept of Logical Domain of Influence (LDOI)41; the LDOI represents the influence of a (set of) fixed node state(s) on the rest of the components in a system. In this section we first present our proposed method to calculate the LDOI for the Boolean threshold functions governing the Campbell et al. model of plant–pollinator community assembly. Then we verify that the simplified logical functions preserve the LDOI and hence can be implemented to further analyze the effect of extinction in plant–pollinator networks. Next, we address one of the main questions that motivated this study: Can stable motif driver set analysis facilitate the identification of keystone species? We discuss the identification of the driver sets of inactive stable motifs and motif groups and present the results of stabilizing these sets to measure the magnitude of the effect of species extinction on the communities. Lastly we discuss possible prevention and mitigation measures based on the knowledge acquired from driver sets of stable motifs and motif groups.Figure 2Illustration of species extinction and restoration in a hypothetical 6-species community. (a) The interaction network (on the left), and the maximal richness community possible for this network (the community with the most established species). Nodes highlighted with green represent established species. (b) The initial extinction of two species, po_1 and po_2 (left) and the community that results after cascading extinctions (right). Nodes highlighted with grey represent extinct species. (c) An intervention to restore pl_2 (left), which induces the restoration of further species, finally leading to a restored community with all the species present (right). The nodes highlighted with teal represent the restored species.Full size imageLDOI in the Boolean threshold modelThe LDOI concept was originally defined on Boolean functions expressed in a disjunctive prime form. Here we extend it to Boolean threshold functions. We implemented it as a breadth first search on the interaction network, as exemplified in Fig. 3. Assume that we want to find the LDOI of a (set of) node(s) (S_0={n_1,dots ,n_N}) and their specific fixed state (Q(S_0)={sigma _{n_1},dots ,sigma _{n_N}}). Starting from the set (S_0), the next set of nodes (S_1) that can acquire a fixed state due to the influence of (Q(S_0)) consists of the nodes that have an incoming edge from the nodes in the set (S_0) in the interaction network. The nodes in set (S_1) are the subject of the first search level. For each node (n_i in S_0) and (n^prime _i in S_1) we assume a “worst case scenario” (i.e., maximal opposition of the effect of (n_i) on (n^prime _i) from other regulators) to find the possible sufficiency relationships between the two. There are five cases:

1.

If (n_i) is a positive regulator of (n^prime _i), then (sigma _{n_i}=1) is a candidate for being sufficient for (sigma _{n^prime _i}=1). We assume that all other positive regulators of (n^prime _i) that have an unknown state (i.e., are not in (Q(S_0))) are inactive and all negative regulators of (n^prime _i) that have an unknown state are active. If (sum _j W_{ij} > 0) under this assumption, then the active state of (n_i) is sufficient to activate (n^prime _i). The virtual node (n^prime _i) that corresponds to (sigma _{n^prime _i}=1) is added to LDOI((Q(S_0))).

2.

If (n_i) is a positive regulator of (n^prime _i), then (sigma _{n_i}=0) is a candidate for being sufficient for (sigma _{n^prime _i}=0). We assume all other positive regulators of (n^prime _i) that have an unknown state are active and all negative regulators of (n^prime _i) that have an unknown state are inactive. If (sum _j W_{ij}le 0) under this assumption, then the inactive state of (n_i) is sufficient to deactivate (n^prime _i). The virtual node (sim n^prime _i) that corresponds to (sigma _{n^prime _i}=0) is added to LDOI((Q(S_0))).

3.

If (n_i) is a negative regulator of (n^prime _i), then (sigma _{n_i}=1) is a candidate for being sufficient for (sigma _{n^prime _i}=0). We assume all positive regulators of (n^prime _i) that have an unknown state are active and all other negative regulators of (n^prime _i) that that have an unknown state are inactive. If (sum _j W_{ij}le 0) under this assumption, then the active state of (n_i) is sufficient to deactivate (n^prime _i). The virtual node (sim n^prime _i) that corresponds to (sigma _{n^prime _i}=0) is added to LDOI((Q(S_0))).

4.

If (n_i) is a negative regulator of (n^prime _i), then (sigma _{n_i}=0) is a candidate for being sufficient for (sigma _{n^prime _i}=1). We assume all positive regulators of (n^prime _i) that have an unknown state are inactive and all other negative regulators of (n^prime _i) that that have an unknown state are active. If (sum _j W_{ij} > 0) under this assumption, then the inactive state of (n_i) is sufficient to activate (n^prime _i). The virtual node (n^prime _i) that corresponds to (sigma _{n^prime _i}=1) is added to the LDOI((Q(S_0))).

5.

If none of the past four sufficiency checks are satisfied, the node (n^prime _i) will be visited again in the next search levels.

The second set of nodes that can be influenced, (S_2), are the nodes that have an incoming edge from the nodes in the set (S_1). The algorithm goes over these nodes in the second search level as described above. This search continues to all the levels of the search algorithm until all nodes are visited (possibly multiple times) and either acquire a fixed state and are added to the LDOI or their state will be left undetermined at the end of the algorithm. In Fig. 3, we illustrate this search to find the LDOI((sim )pl_1). The first search level is (S_1={)po_1, po_3(}); (sim )pl_1 is sufficient to deactivate po_3, but not po_1. As a result, (sim )po_3(in ) LDOI((sim )pl_1). This process continues until all levels are visited and at the end of the algorithm LDOI((sim )pl_1()={sim )po_3, (sim )pl_2, (sim )pl_3, (sim )pl_4, (sim )pl_5, (sim )po_1, (sim )po_2 (}).Figure 3Breadth first search of the interaction network to find the LDOI of a (set of) fixed note state(s) in Boolean threshold functions governing the dynamics of plant–pollinator networks. (a) An interaction network with five plants and 3 pollinators. (b) The breadth first search in the case of starting from the node state (sim )pl_1. The nodes with incoming edges from pl_1 make up (S_1={)po_1, po_3(}). The second sufficiency check is satisfied for node state (sim )po_3, as a result (sim )po_3(in ) LDOI((sim )pl_1). The same process is applied for node po_1, but none of the sufficiency checks are satisfied, so this node will be visited again later. The next level of the search consists of the nodes that have incident edges from (S_1), i.e., (S_2={)pl_2, pl_3, pl_4, pl_5(}). The second sufficiency check is satisfied for all of these nodes and they are all fixed to their inactive state in the LDOI((sim )pl_1). Lastly, we reach (S_3={)po_1, po_2(}). Node po_1 is reached again, and with both its positive regulators fixed to their inactive states the second sufficiency check is satisfied and node po_1 is fixed to its inactive state as well. The same holds for po_2 and hence LDOI((sim )pl_1()={sim )po_3, (sim )pl_2, (sim )pl_3, (sim )pl_4, (sim )pl_5, (sim )po_1, (sim )po_2 (}).Full size imageTo measure the accuracy of the simplification method originally introduced in28, we analyzed logical domains of influence in 6000 networks with 50–70 nodes. These networks are among the largest in our ensembles and have the most complex structures. We randomly selected (sets of) inactive node states, found their LDOIs using the Boolean threshold functions and the simplified Boolean functions, and compared the two resulting LDOIs. We used 8 single node states and 8 combinations of size 2 to 4 for each network. We found that in all cases the LDOI calculated using the simplified Boolean functions matches the LDOI calculated using the Boolean threshold functions.Next, we analyzed (sets of) active node states and their LDOIs in the same ensembles of networks. Similar to the previous analysis, we used 8 single node states and 8 combinations of size 2 to 4 for each network. Our analysis shows that in 77.1% of the cases the LDOI calculated using the simplified Boolean functions matches the LDOI calculated using the Boolean threshold functions. In 22% of the cases the LDOI calculated from the simplified Boolean functions contains the LDOI calculated from the threshold functions, and it also contains extra active node states, overestimating the LDOI by 57.5% on average. These additional members of the LDOI result from the fact that the simplified Boolean functions contain fewer negative regulators than the threshold functions. The guiding principle of the simplification method is that the probability of (H(x)=1) conserves the probability of each node having an active state across all the states it can have. In contrast, the probability of the propagation of the active state is not necessarily preserved and tends to be higher in the simplified Boolean model; thus the LDOI of the active node states is overestimated in some cases.In the rest of the cases (about 1%), the LDOI calculated from the simplified Boolean functions does not fully capture the LDOI calculated from the threshold functions. This again is caused by the sparsification of the negative edges in the simplified Boolean functions. In the threshold functions, the activation of 4 or more negative regulators of a target node combined with one active positive regulator is sufficient to deactivate the target node, i.e., there might be inactive node states in the LDOI of a set of active node states. However, some of these negative regulators drop in the simplified Boolean model and the inactive state of the target node is not necessarily in the LDOI of the set of active node states in the simplified case. This is the rare mechanism by which the simplified model might underestimate the influence of active node states on the rest of the network.In the following section we are interested in analyzing the effect of species extinction on the established community, i.e., we look at the LDOI of (sets of) inactive node states. Observing that the influence of extinction of species is measured correctly in the simplified Boolean models, we conclude that these models can be utilized to further analyze the process of extinction and its ecological implications.Stable motif based identification of species whose loss leads to cascading extinctionsEach stable motif or motif group can have multiple driver sets; stabilization of each driver set leads to the stabilization of the whole motif or motif group. In plant–pollinator interaction networks, the stable motifs either represent a sub-community (when the constituent nodes stabilize in their active states) or the simultaneous extinction of all species in the group (when the constituent nodes stabilize to their inactive states). Stabilization of the nodes in the driver set of an inactive stable motif results in stabilization of all the nodes in the stable motif to their inactive state, i.e., cascading extinction of the constituent species.The knowledge gained from stable motif analysis and the network of functional relationships offers insight into the cascading effect of an extinction that constitutes a driver set of an inactive stable motif. The magnitude of this effect depends on (i) the number of nodes that the inactive stable motif contains and (ii) the number of virtual nodes (including motifs and motif groups) corresponding to inactive species that are logically determined by the stabilization of the inactive stable motif.To investigate the role of stable motifs in the study of species extinction in plant–pollinator networks, we simulated extinctions that drive inactive stable motifs in 6000 networks with the sizes of 50–70 nodes. We considered driver sets of size 1, 2, or 3, and implemented them by fixing the corresponding node(s) to its (their) inactive state. As a point of comparison, we also performed a “control” analysis using the same networks with the same size of initial extinction; however, the candidates of initial extinction are inactive node states that do not drive stable motifs or motif groups. Based on the properties of the drivers of stable motifs, one expects that following the extinction of driver species, cascading extinctions of other species follow, while the same does not necessarily hold for non-driver species. As a result, we expect to observe greater damage to the original community when driver species become extinct.We assume that the “maximal richness community”—the community (attractor) in which the largest number of species managed to establish—is the subject of species extinction. This maximal richness community results from the stabilization of all active stable motifs. All other attractors that have some established species contain a subset of all active stable motifs and thus will contain a subset of the species of the maximal richness community. While for a generic Boolean model with multiple attractors one expects that a perturbed version of the model also has multiple attractors, this specific perturbation of a plant–pollinator model (namely, extinction of species in the maximal richness community) has a single attractor. We prove this by contradiction. Assume there are two separate attractors in the perturbed model, which means that there is at least one node that has opposite states in these two attractors. Note that this bi-stability is the result of the perturbation and not a property of the original system as the maximal richness community (an attractor) is the starting point for the introduced extinction. Specifically, the inactive state of the extinct node has to lead to the stabilization of another node to its active versus inactive states in the two separate attractors. The only case in which the stabilization of an inactive node state can result in the stabilization of an active node state is if there is a negative edge from the former to the latter in the interaction network after simplification. Since the Boolean function in 2 is inhibitor dominant, the negative regulators that remain in the Boolean model must be in their inactive states in the maximal richness attractor. As they are already inactive (extinct), they are not candidates for extinction. The only nodes that are candidates for extinction are the ones that positively regulate other nodes; perturbing the system by fixing these candidates to their inactive states cannot lead to the active state of a target node. In conclusion, bi-stability is not possible.We found the new attractor of the system given the (combination of) inactive node state(s) using the the functions percolate_and_remove_constants() and trap_spaces() from the pyboolnet Python package. We quantify the effects of the initial extinction(s) on the maximal richness attractor by the percentage change in the number of active species, which we call damage percentage. Note that this choice of maximal richness community as the reference and starting point allows us to detect the cascading extinctions following the initial damage.In Fig. 4 the left column plots show the average damage percentage caused by the extinction of 1 (top panel), 2 (middle panel), or 3 (bottom panel) species that represent driver sets of stable motifs and motif groups, while the right column plots illustrate the average damage percentage as a result of the extinction of 1, 2 or 3 species that represent non-driver nodes. Comparing the two columns, one can notice that stabilization of the driver sets of stable motifs and motif groups leads to considerably larger damage to the communities. This is due to the fact that stabilization of driver sets ensures the stabilization of entire inactive stable motifs and motif groups and hence ensures cascading extinctions. Comparing the plots in the left column, we see that the larger the driver sets are, the larger the damage to the community becomes. This is because larger driver sets are more likely to stabilize larger stable motifs and motif groups. This figure illustrates the significance of stable motifs and their driver sets in the study of species extinction in plant–pollinator communities.Figure 4Histogram plots illustrating the average percentage of the damage caused in an established community after the extinction of species. This analysis is performed over 6000 networks with the size of 50–70 nodes. To study the extinction of species we started from the maximal richness community, then we fixed the nodes that correspond to the focal species to the their inactive states. The original extinctions are excluded from the damage percentages. The left column plots show the average damage percentage caused to the maximal richness community by the extinction of a driver set of size 1 (top), 2 (middle), or 3 (bottom) of an inactive stable motif or motif group. For each network, we determined all the relevant driver sets of one stable motif or motif group, we performed the extinction and calculated the resulting damage, then we calculated the average damage percentage over all data points collected for the same network. The right column plots show the average damage percentage caused to the maximal richness community by the extinction of 1 (top), 2 (middle), and 3 (bottom) non-driver, randomly chosen nodes. Each time a randomly selected combination of non-driver nodes were the subject of simultaneous extinction until all combinations are explored and then we calculated the average damage percentage over all data points collected for each network. The number of networks that qualify for each of these 6 categories differ (e.g., some networks have a stable motif with a driver set of size 2 but no stable motif with a driver set of size 3). In the left column 5529, 3212, and 1980 networks and in the right column 5779, 5626, and 5423 networks qualified respectively. The red lines represent the mean value of all the presented data points in each plot.Full size imageIn Fig. 4 left column, the full driver set of one inactive stable motif or motif group was stabilized. However, the species that become extinct might only contain a subset of a driver set of a stable motif or motif group, i.e., they only stabilize a subset of the inactive node states in the stable motif or motif group. We compare the extinction effect caused by the stabilization of a full driver set of four nodes with the effect of the extinction of four nodes that contain a partial driver set in Fig. 5 using the batch of the largest networks in this study, i.e, the batch that contains networks with 30 nodes representing plant species and 40 nodes representing pollinator species. This choice is due to the fact that the existence of stable motifs and motif groups having a driver set of four node states is highly probable in larger networks. As expected, the stabilization of the complete driver set leads to greater damage. Stabilization of the same number of nodes that contain a partial driver set leads to significantly less damage and species loss in the community; the median damage percentage in the case of stabilization of partial driver sets is 22.6% while it is 69.2% in the case of stabilization of the full driver sets. We also note that damage of more than 90% occurs rarely and is only possible when a full driver set is stabilized (see Fig. 5 right plot). This suggests that the motif groups that lead to total extinction tend to have a driver set with more than four nodes; in other words, only the simultaneous extinction of five or more species would lead to total community collapse.Figure 5Histogram plots illustrating the average percentage of the damage caused in an established community after the extinction of species. This analysis is performed over 1000 networks with the size of 70 nodes (30 nodes representing plant species and 40 nodes representing pollinator species). The original extinctions are excluded from the damage percentages. The left plot shows the average damage percentage caused to the maximal richness community by the extinction of 2 species that are a subset of the 4-node driver set of an inactive stable motif or motif group plus 2 randomly selected non-driver species. The right plot shows the damage percentage caused to the maximal richness community by the extinction of 4-node driver sets of the same inactive stable motifs and motif groups. Each time the driver set of one stable motif or motif group was the subject of extinction and we calculated the average damage percentage over all data points collected for each network. 295 networks qualified for this analysis.Full size imageMotif driver set analysis outperforms structural measures in identifying keystone speciesThe literature on ecological networks offers multiple measures that reflect the importance of each species for community stability. One family of such measures is centrality (quantified by the network measures degree centrality and betweenness centrality). Previous studies45,46 have shown that species (nodes) with higher centrality scores are keystone species in ecological communities (i.e., species whose loss would dramatically change or even destroy the community). The nodes with highest in-degree centrality (such as pl_2 in Fig. 6a) represent generalist species that can receive beneficial interactions from multiple sources and survive. The nodes with highest betweenness centrality (such as pl_2 and po_2 in Fig. 6a) represent species that act as connectors and help the community survive. We find that high centrality corresponds to specific patterns in the expanded network: the inactive state of generalist or connector species is often the driver of a cascading extinction. Indeed, stable motif analysis of the expanded network in Fig. 6b confirms that there is an inactive stable motif (highlighted with grey) driven by the minimal set {(sim )pl_2}. The fact that node pl_2 is a stable motif driver means that in the case of the extinction of pl_2 the whole community collapses.To compare the effectiveness of stable motif analysis to the effectiveness of the more studied structural measures to identify keystone species, we performed an analysis similar to the previous section. We compared the magnitude of cascading extinctions in the case of extinction of stable motif driver nodes and of nodes with high values of previously introduced structural importance measures. Specifically, we used node betweenness centrality, node contribution to nestedness47, and mutualistic species rank (MusRank)22 to find crucial species based on their structural properties. For more details on definition and adaptation of these two measures see “Methods”. In this analysis, we used each measure to target species in the simplified Boolean models as follows:

1.

Betweenness centrality: The 10% of species with the highest betweenness centrality are chosen to be candidates for extinction.

2.

Node contribution to nestedness: The species with the most interactions tend to contribute the least to the community nestedness. Targeting them most likely leads to a faster community collapse48. As a result, 10% of species with the lowest contribution to network nestedness are chosen to be candidates for extinction. For more details on this measure, please see “Methods”.

3.

Pollinator MusRank: The pollinator species with the highest MusRank importance are more likely to interact with multiple plants, so the 10% of pollinator species with the highest importance are chosen to be candidates for extinction. For more details on this measure, please see “Methods”.

4.

Plant MusRank: The plant species with the highest MusRank importance are more likely to interact with multiple pollinators, so the 10% of plant species with the highest importance are chosen to be candidates for extinction.

Figure 7 illustrates the results of this analysis in 6000 networks with 50–70 nodes. In each network the 1-node, 2-node, and 3-node driver sets of inactive stable motifs are identified and made extinct. In the same networks 10% of nodes based on betweenness centrality, node contribution to nestedness, and node MusRank score were chosen to be candidates for extinction. To match the “driver set” data, all choices of 1, 2, or 3 nodes in these sets were explored and the damage was averaged over each extinction size for each network. We observe the cascading extinction and calculate the damage percentage relative to the maximal richness attractor. The plot represents the collective data over all initial simultaneous extinction sizes of 1, 2, and 3 species.Comparing the four methods, one notices that the histograms acquired using stable motif driver sets, node betweenness centrality, and node contribution to nestedness are very similar, showing a peak for the 10–20% bin of the damage, and a long tail that reaches a damage percentage of 80–100%. The MusRank score performs less well in identifying the crucial species. Also, the frequency of the higher damage percentages shows that node contribution to nestedness is the closest to the “driver set” method in identifying nodes whose extinction causes the collapse of the whole community, making it the best structural measure out of the three. Nevertheless, the driver set method finds keystone species that cannot be identified via structural measures, as the corresponding damage percentage histogram has the most prominent tail at the right edge of the panel. Indeed, stable motif driver sets identified 82%, 80%, and 546% more species whose extinction leads to 60% or higher damage to the community when compared to betweenness centrality, node nestedness, and node MusRank score based methods respectively.The reason for the higher effectiveness of driver set analysis is illustrated in Fig. 8 in which the MusRank score and node contribution to nestedness are calculated for two example networks. One can see how these two measures might incorrectly identify less vital species. In the left column of Fig. 8, MusRank identifies the node po_2, highlighted with green, as the most important species. However, this node does not have any outgoing edges; its extinction does not lead to any cascading extinction. The inability of the MusRank score to consider the direction of edges causes such misidentification. In the right column, the three nodes highlighted with yellow have the lowest contributions to network nestedness. The expanded network shows that these three nodes together are not able to cause full community collapse, while the three-node driver set of the inactive stable motif can. Since the nestedness definition depends on the number of mutual interactions, it might fail to identify some of the keystone nodes that are necessary to the stability of the community (for more details on node nestedness see “Methods”).Previously it was shown that identifying the stable motifs and their driver sets can successfully steer the system toward a desired attractor or away from unwanted ones37,38,43. Stable motif analysis of the Boolean model offers insight into the dynamical trajectories of the system; hence control strategies can be developed accordingly. In the next section we use stable motif driver sets to suggest control methods and analyze their efficiency.Figure 6Generalist species in the interaction network and the expanded network. (a) A simplified network consisting of 3 plant and 3 pollinator species. pl_2 is a generalist species, i.e., it has two incoming edges indicating that it can survive on either of its sources of pollination, po_1 or po_2. The expanded network in (b) illustrates that the stabilization of the grey stable motif stabilizes all the nodes to their inactive states, and hence causes full community collapse. (sim )pl_2 is the minimal driver set of the grey stable motif, consistent with the strong damage induced by the loss of a generalist species.Full size imageFigure 7Histogram plots illustrating the performance of driver set analysis versus structural measures in identifying keystone species. The analysis was done on 6000 networks with sizes of 50–70 nodes. The starting point is the maximal richness community, i.e., the attractor in which the most species establish. For each network 1, 2, and 3 node(s) were selected and simultaneously fixed to their inactive states. After the cascading damage the new attractor is compared to the maximal richness attractor to calculate the damage percentage. The structural measures—betweenness centrality, node nestedness contribution, and node MusRank score—were calculated for all nodes in each network; the top 10% according to the relevant ordering were candidates to being fixed to their inactive states. The network IDs were matched, i.e., only the networks that had candidate nodes according to all four measures for each extinction size are included in this plot. The total number of data points is 6360. The red solid lines represent the mean and the black dashed lines represent the median over all data points in each plot.Full size imageFigure 8Networks illustrating examples of when structural measures fail to identify keystone species. In both columns simplified networks consisting of 3 plant and 3 pollinator species are presented. The MusRank is calculated for all the nodes in the network in the left column and denoted in the node labels. The expanded network corresponding to this network is shown below. Node contribution to network nestedness is calculated for all the nodes in the network in the right column and denoted in the node labels. Similarly the expanded network that correspond to it is shown below. Note that these two networks have different edges. In the left column MusRank score identifies node po_2, highlighted with green, as the most important, while the expanded network shows that the extinction of po_2 does not cause any further damage to the community, as this node has no outgoing edges. This is due to the fact that MusRank calculation process fails to consider the directed network and replaces all the directed edges with undirected ones. The MusRank score does not identify po_3 as a crucial species; however, virtual node (sim )po_3, outlined with black in the expanded network is a driver of a stable motif that has all other nodes in its LDOI; the extinction of po_3 leads to full community collapse. In the right column, the nodes highlighted with yellow (pl_2, pl_3, and po_2) have the lowest node contribution to nestedness, which predicts that these nodes are likely crucial to the stability of the community. Analyzing the expanded network, one can see that these three nodes together are not able to drive the inactive stable motif highlighted with teal. The minimal driver set for this stable motif, outlined with black, consists of {(sim )po_1, (sim )po_2, (sim )po_3}; together these nodes drive the inactive stable motif and cause full community collapse. The nestedness-based measure was not able to capture the significance of nodes po2 and po_3.Full size imageDamage mitigation measures and strategies for endangered communitiesThere are two substantial questions related to managing the damage induced by species extinction: (1) How can one prevent the damage as much as possible? (2) Once the damage happens, the reintroduction of which species can restore the community and to what extent? In this section we aim to answer these questions in the context of the Campbell et al. model, implementing stable motif based network control. This analysis can inform agricultural and ecological strategies employed to prevent and mitigate damage.Damage preventionOne of the most important questions in ecology is what strategies to use so that we can prevent and avert extinction damage to the community. In this section we analyze how the knowledge from stable motif analysis and driver sets can be implemented to minimize the effect of extinction of keystone species in case of limited resources. Each attractor of the original system can have multiple control sets; stabilizing the node states in each control set ensures that the system reaches that specific attractor. The same information from the attractor control sets can be implemented to prevent the system from converging into unwanted attractors. Zañudo et al. illustrated that by blocking (not allowing to stabilize) the stable motifs that lead to the unwanted attractors, one can decrease the probability (sometimes to zero) that the system arrives in those attractors38. In order to block an attractor, the control sets of that attractor are identified and the negations of the node states in the control sets are externally imposed. This approach eliminates the undesired attractor; however, new attractors might form that are similar to the eliminated attractor. Campbell et al. showed that in order to avoid such new attractors one needs to block the parent motif, which in this case is the largest strongly connected subgraph of the expanded network that contains the inactive virtual nodes44. Here, we investigate how stable motif blocking based attractor control can identify the species whose preservation would offer the highest benefit in avoiding catastrophic damage to the community. This information would aid the development of management strategies in plant–pollinator communities.To avoid all attractors that lead to some degree of species extinction, one needs to block all the driver sets of all inactive stable motifs and motif groups in a given network. Implementing this in 100 randomly selected networks with 25 plant and 25 pollinator nodes, we found that 45.6% of the species in the maximal richness community need to be kept (prevented from extinction) to ensure the lack of cascading extinctions. Given that management resources are usually limited, active monitoring and conservation of almost half of the species in a community seems costly and impractical. Hence, we set a more feasible goal of identifying and blocking the driver set(s) of the largest inactive stable motif or motif group in each network. The same 100 networks containing 50 nodes are the subject of analysis in this section. The reason for performing the analysis in a relatively limited ensemble is that it involves the identification of all driver sets of the largest inactive stable motif or motif group, which is computationally expensive. For each network, the driver set of the largest inactive stable motif or motif group (which corresponds to the extinction of all the species in that group) is identified and blocked (that is, the corresponding species are not allowed to go extinct). Then the same number of species as in the driver set of that stable motif or motif group are selected and stabilized to their inactive state. We considered all combinations of node extinctions outside the blocked subset, calculated the damage percentage relative to the maximal richness community, and then averaged over all data points for each network. As a control, we repeated the analysis without blocking; the size of the initial extinction is the same as in the previous analysis for consistency.Figure 9 shows the result of the analysis described above for 100 networks. The left box and whiskers plot illustrates the damage percentage relative to the maximal richness community when the blocking feature is activated, while the right box and whiskers plot shows the damage percentage relative to the maximal richness community when the blocking is disabled. The average and median damage percentages are 14.96% and 13.04% respectively when the largest inactive stable motif or motif group was blocked and 24.73% and 20.38% when it was not. This (sim )10% difference in the average between the two sets of results, as well as the fewer cases of high-damage outliers in the left plot, demonstrates that by preventing the extinction of species identified by stable motif analysis, one can prevent catastrophic community damage considerably.To estimate the fraction of species that would need to be monitored to prevent their extinction, we compared the size of the maximal richness attractor and the size of the driver set of the largest stable motif. The maximal richness community represents an average of 32% of the original species pool, approximately 15 out of 50 species. The driver sets of the largest stable motifs had an average size of 2.5 node states over all 100 networks, i.e., about 16.6% of the maximal richness community. In ecological terms, given limited resources, the information gained from stable motif driver sets can help direct the conservation efforts toward the keystone species that play a key role in maintaining the rest of the community in a cost-effective manner.Figure 9Box plots comparing the damage communities face if the largest inactive stable motif or motif group is completely blocked, i.e., all the drivers of this inactive stable motif or motif group are prevented from stabilizing versus if the same stable motif or motif group is allowed to stabilize. This analysis was performed over 100 randomly selected networks that contain 25 plant and 25 pollinator nodes. All the driver sets of an inactive stable motif or motif group are identified. From left to right the box and whiskers plots show the average damage percentage relative to the maximal richness community if the largest inactive stable motif is blocked and the same quantity if the largest stable motif or motif group is not blocked respectively. For the left box and whiskers plot, all combinations of inactive node states except the driver sets are considered, and for the right box and whiskers plot all combinations are explored. Due to the computational complexity caused by combinatorial explosion, this analysis was performed over 100 randomly selected 50-node networks.Full size imageRestoration of a group of speciesAlthough human preservation efforts have been directed toward community conservation, there are many industrial activities that lead to ecosystem degradation. Ecologists are interested in developing restoration strategies to be deployed after a stable community is hit by catastrophic damage to recover biodiversity and the ecosystem functions it provides49. Here we propose that stable motif analysis and the driver sets identified from the expanded network can give insight into restoration measures. While we examined the inactive stable motifs in the study of species extinction, here we focus on the active stable motifs as our goal is to restore as much biodiversity as possible.Several network measures have been proposed to identify the species that if re-introduced would restore the community considerably. Two of the most studied algorithms include maximising functional complementarity (or diversity) and maximising functional redundancy50. The first strategy targets the restoration of the species that provide as many functions to the ecosystem as possible; this approach results in a community that has a maximal number of functions provided by different groups of species. Alternatively, maximising the functional redundancy yields a community in which several species perform the same function. While this resultant community might have a limited number of functions, it is robust. Both of these community restoration approaches have been studied extensively (e.g. see21).We hypothesize that restoring the species that constitute driver sets of active stable motifs can help maximise the number of species post-restoration. Since there is evidence that functional diversity correlates with the number of species in the community51, we compare the post-restoration communities identified by stable motif driving with the functional diversity maximisation approach. As discussed in section LDOI in the Boolean threshold model, the Boolean simplification of the threshold functions leads to an overestimation of the LDOI of active node states (compared to the original threshold functions) in some networks. We evaluate the negative effects of this overestimation by checking the effectiveness of the restored species in the original threshold model.The same 6000 networks we examined in the last section were the subject of this analysis. To create an unbiased initial community, we create the damaged communities by eliminating the same number of species from the maximal richness community as the number that will be restored. We identify the inactive stable motif or motif group with the driver set size of 1, 2, or 3 node states that causes the most damage to the maximal richness community. We then eliminate the species corresponding to this driver set to reach the most damaged community for the given size of the initial extinction. This community is the starting point for two analyses. In the stable motif driving approach we stabilized an active stable motif that has a driver set of the same size as the initial extinction to reach a post-restoration community and calculated the percentage of the extinct species that were restored. In the functional diversity maximization based approach we re-introduced the same number of species selected from the to 10% of species in terms of their contribution to functional diversity.To calculate the functional diversity of a community one needs to (1) define and construct a trait matrix, (2) determine the distance (trait dissimilarity) of pairs of species, (3) perform hierarchical clustering based on the distances to create a dendrogram, and (4) calculate the total branch length of the dendrogram, i.e., the sum of the length of all paths51,52. Petchey et al. argued that resource-use traits among plant and pollinator species can be used to classify the organisms into separate functional groups53 and Devoto et al. proposed the use of the adjacency matrix based on the interaction network as the trait matrix21. In this study we do the same and implement the bipartite adjacency matrix to construct the distance matrix.Since the networks of the Campbell et al. model are directed, we modify the algorithm in that we have two separate adjacency matrices, one denoting the edges incoming to plant species and the other denoting the edges incoming to pollinator species. The hierarchical clustering algorithm is then run on each of these matrices separately, resulting in a dendrogram for each adjacency matrix. If extinction occurs in a community, the functional diversity of the survived community can be determined by calculating the total branch length of the subset of the dendrogram that includes only the survived species. The restoration strategy using this method is to re-introduce the nodes whose branches add the most to the total branch length of this subset, i.e., maximise the functional diversity of the survived community54. For more details see “Methods”.In each network, the percentage of the extinct species that were restored was calculated and averaged over all data points for each restoration size and each network. Figure 10 illustrates the results of this investigation. Applied to the simplified Boolean model, the median restoration percentage in the case of active stable motif driver set method (blue plot) is 80%. The functional diversity maximization strategy to restoration (yellow plot) yields a lower median restoration percentage, 73%, as well as a large number of low-restoration outliers. Although one might argue that identifying beneficial species using the functional diversity maximization strategy works well, the higher percentage of the cases of 80–100% restoration in case of the active stable motif driver set analysis indicates that the latter identifies some of the most effective restorative species that are not identified via the former method. As in a minority of cases the simplified Boolean model overestimates the positive impact of the sustained presence of a species (see section LDOI in the Boolean threshold model), we sought to verify the effectiveness of the predicted restoration candidates in the original threshold model. The blue (respectively, yellow) box and whiskers plot on the right represents the restoration percentages of the same species as in the left blue (respectively, yellow) plot when these species are restored in the threshold model. The median of the right blue plot is 70%, while the median of the right yellow is 63%, preserving the advantage of the stable motif driver sets. We conclude that although the simplified Boolean model overestimates the restoration effectiveness of certain driver sets (visible in the fact that the lower whisker of the blue plot on the right goes well below the lower whisker of the blue plot on the left), stable motif driver sets are more effective in both comparisons.Figure 10Box and whiskers plots illustrating the average percentage of the extinct species that are restored following the stable motif driver set restoration strategy (blue) versus the functional diversity based approach (yellow). This analysis is performed over 6000 networks with sizes of 50–70 nodes. Starting from the maximal richness community, for each network one inactive stable motif with a driver set of 1, 2 or 3 nodes was stabilized to reach a new damaged community. This task was performed until the community with the most extinct species was identified. This is the community we set as the starting point for the restoration process using both methods. The pair on the left represents the two methods applied to the simplified Boolean model. For both methods we identified 1, 2, or 3 influential nodes for community restoration and we calculated the percentage of the extinct species that could be restored. The pair on the right represents restoring the same species identified by each method in the previous analysis in the original threshold model. In all analyses the community restoration percentage was averaged over all combinations of the same size, for each network and each method. The IDs of all networks are matched.Full size imageCommunity restoration via attractor controlAs illustrated in section “Restoration of a group of species”, stable motif analysis identifies promising and cost-effective group restoration strategies. In this section we aim to go further and identify interventions that can maximally restore a community. Previous stable motif based network control methods37,38,55 require a search for the smallest set of node states to control the system once the stable motif stabilization trajectories are identified. This smallest set may not contain a node from each stable motif in the sequence. In this work, however, we know that each stable motif or motif group needs to be controlled individually28 because the stabilization of none of the motifs results in the stabilization of another. As a result, the control set of each attractor is the same as the union of the driver sets of all members in the consistent combination corresponding to that attractor.In this section we examined this attractor control method by setting the communities with 70% or more of the species in the maximal richness community as the target, i.e., the attractors that have 70% of the species in the maximal richness community are assumed to be the desired attractors. We then recorded the size of the minimal control set needed to achieve each of these attractors. Note that stabilizing each of these control sets guarantees that the system reaches the corresponding attractor38.For this section, we analyzed 6000 networks that have 50–70 nodes. Figure 11 represents box-and-whiskers plots of the size of the minimal set of species that need to be restored, where the target community sizes are classified into three groups based on the percentage of the species relative to the maximal richness attractor. One can see that in half of the cases, the restoration of either 1 or 2 species manages to restore more than 70% of the maximal richness community. The largest set has 8 species that need to be restored; however, this data point is an outlier. As illustrated, driver set analysis and stable motif based attractor control can efficiently identify the species that play an influential restorative role and suggest management strategies that are effective at the scale of the whole community. To assess the impact of the LDOI inflation on this result, we used the restoration candidates identified by control sets of the attractors of the Boolean model in the threshold functions of a subset of networks. The results of comparing the restoration percentage is shown in Fig. 14. The first quartile, median and third quartile values are 78.26%, 86.6%, and 100% for the simplified Boolean models and 43.78%, 72.41%, and 85.71% for the threshold model.To further compare the results of restoration obtained from the two models we sorted the species in the order of their contribution to community restoration following a catastrophic damage. We randomly selected 100 of the largest (70-node) networks, which have the highest probability of a discrepancy between the threshold functions and the simplified Boolean model. In 72% of the cases the two rankings matched completely, and in the majority of the remaining cases only one species was misplaced in the simplified Boolean model-based ranking. To conclude, there is a significant advantage to the implementation of the simplified Boolean model and the drawback can be addressed by a follow-up checking on the original threshold functions.Figure 11The number of species that need to be restored to save 70% of more of the species in the maximal richness community. In this analysis 6000 networks with 50–70 nodes were the subject. For each networks all the attractors that have 70% or more of the species in the maximal richness attractor are identified and set to be the target attractors. The control set of these attractors are then classified into three groups based on the percentage as illustrated in the figure. From left to right, the box and whiskers represent the size of the control set of attractors that have 70–80%, 80–90%, and 90–100% of the species in the maximal richness attractor respectively.Full size image More

Plaisance, L., Caley, M. J., Brainard, R. E. & Knowlton, N. The diversity of coral reefs: what are we missing? PLoS One. 6, e25026 (2011).Article 
ADS 
CAS 

Google Scholar 
Hoegh-Guldberg, O., Poloczanska, E. S., Skirving, W. & Dove, S. Coral reef ecosystems under climate change and ocean acidification. Front. Mar. Sci. 4, 158 (2017).Article 

Google Scholar 
Mora, C. et al. Global human footprint on the linkage between biodiversity and ecosystem functioning in reef fishes. PLoS Biol. 9, e1000606 (2011).Article 
CAS 

Google Scholar 
Burke, L., Reytar, K., Spalding, M. & Perry, A. Reefs at Risk Revisited. 130 pp. (World Resources Institute, Washington, D.C., 2011).Hicks, C. C., Graham, N. A. J., Maire, E. & Robinson, J. P. W. Secure local aquatic food systems in the face of declining coral reefs. One Earth. 4, 1214–1216 (2021).Article 
ADS 

Google Scholar 
Cinner, J. E. et al. Gravity of human impacts mediates coral reef conservation gains. PNAS 115, E6116–E6125 (2018).Article 
CAS 

Google Scholar 
Eddy, T. D. et al. Global decline in capacity of coral reefs to provide ecosystem services. One Earth. 4, 1278–1285 (2021).Article 
ADS 

Google Scholar 
Graham, N. A. J. et al. Human disruption of coral reef trophic structure. Curr. Biol. 27, 231–236 (2017).Article 
CAS 

Google Scholar 
Sherman, C. S., Heupel, M. R., Moore, S. K., Chin, A. & Simpfendorfer, C. A. When sharks are away rays will play: effects of top predator removal in coral reef ecosystems. Mar. Ecol. Prog. Ser. 641, 145–157 (2020).Article 
ADS 

Google Scholar 
Ruppert, J. L. W., Travers, M. J., Smith, L. L., Fortin, M.-J. & Meekan, M. G. Caught in the middle: combined Impacts of shark removal and coral loss on the fish communities of coral reefs. PLoS One. 8, e74648 (2013).Article 
ADS 
CAS 

Google Scholar 
Last, P. R. et al. Rays of the World. (CSIRO Publishing, 2016).Ebert, D. A., Dando, M. & Fowler, S. Sharks of the World. 2nd edn, 608 (Princeton University Press, 2021).Heupel, M. R., Lédée, E. J. I. & Simpfendorer, C. A. Telemetry reveals spatial separation of co-occurring reef sharks. Mar. Ecol. Prog. Ser. 589, 179–192 (2018).Article 
ADS 

Google Scholar 
Heupel, M. R., Papastamatiou, Y. P., Espinoza, M., Green, M. E. & Simpfendorfer, C. A. Reef shark science – key questions and future directions. Front. Mar. Sci. 6, 12 (2019).Article 

Google Scholar 
Roff, G., Brown, C. J., Priest, M. A. & Mumby, P. J. Decline of coastal apex shark populations over the past half century. Commun. Biol. 1, 223 (2018).Article 

Google Scholar 
Williams, J. J., Papastamatiou, Y. P., Caselle, J. E., Bradley, D. & Jacoby, D. M. P. Mobile marine predators: an understudied source of nutrients to coral reefs in an unfished atoll. Proc. R. Soc. B. 285, 20172456 (2018).Article 

Google Scholar 
Heithaus, M. R., Wirsing, A. J. & Dill, L. M. The ecological importance of intact top-predator populations: a synthesis of 15 years of research in a seagrass ecosystem. Mar. Freshw. Res. 63, 1039–1050 (2012).Article 

Google Scholar 
Peel, L. R. et al. Stable isotope analyses reveal unique trophic role of reef manta rays (Mobula alfredi) at a remote coral reef. R. Soc. Open Sci. 6, 190599 (2019).Article 
ADS 
CAS 

Google Scholar 
O’Shea, O. R., Thums, M., van Keulen, M. & Meekan, M. Bioturbation by stingrays at Ningaloo Reef, Western Australia. Mar. Freshw. Res. 63, 189–197 (2012).Article 

Google Scholar 
Takeuchi, S. & Tamaki, A. Assessment of benthic disturbance associated with stingray foraging for ghost shrimp by aerial survey over an intertidal sandflat. Continental Shelf Res. 84, 139–157 (2014).Article 
ADS 

Google Scholar 
Burkholder, D. A., Heithaus, M. R., Fourqurean, J. W., Wirsing, A. & Dill, L. M. Patterns of top-down control in a seagrass ecosystem: could a roving apex predator induce a behaviour-mediated trophic cascade? J. Anim. Ecol. 82, 1192–1202 (2013).Article 

Google Scholar 
Creel, S. & Christianson, D. Relationships between direct predation and risk effects. TRENDS Ecol. Evolution. 23, 194–201 (2008).Article 

Google Scholar 
Ward-Paige, C. A. et al. Large-scale absence of sharks on reefs in the greater-Caribbean: a footprint of human presence. PLoS One. 5, e11968 (2010).Article 
ADS 

Google Scholar 
Espinoza, M., Cappo, M., Heupel, M. R., Tobin, A. J. & Simpfendorfer, C. A. Quantifying shark distribution patterns and species-habitat associations: implications of marine park zoning. PLoS One. 9, e106885 (2014).Article 
ADS 

Google Scholar 
Graham, N. A., Spalding, M. D. & Sheppard, C. R. Reef shark declines in remote atolls highlight the need for multi-faceted conservation action. Aquat. Conserv.: Mar. Freshw. Ecosyst. 20, 543–548 (2010).Article 

Google Scholar 
Nadon, M. O. et al. Re-creating missing population baselines for Pacific reef sharks. Conserv. Biol. 26, 493–503 (2012).Article 

Google Scholar 
MacNeil, M. A. et al. Global status and conservation potential of reef sharks. Nature 583, 801–806 (2020).Article 
ADS 
CAS 

Google Scholar 
Dulvy, N. K. et al. Overfishing drives over one-third of all sharks and rays toward a global extinction crisis. Curr. Biol. 31, 1–15 (2021).Article 

Google Scholar 
Walls, R. H. L. & Dulvy, N. K. Eliminating the dark matter of data deficiency by predicting the conservation status of Northeast Atlantic and Mediterranean Sea sharks and rays. Biol. Conserv. 246, 108459 (2020).Article 

Google Scholar 
Yan, H. F. et al. Overfishing and habitat loss drives range contraction of iconic marine fishes to near extinction. Science Adv. 7, eabb6026, (2021).Butchart, S. H. M. et al. Using Red List Indices to measure progress towards the 2010 target and beyond. Philos. Trans. R. Soc. B 360, 255–268 (2005).Article 
CAS 

Google Scholar 
Sherman, C. S. et al. Taeniura lymma. The IUCN Red List of Threatened Species, eT116850766A116851089 (2021). 10.2305/IUCN.UK.2021-1.RLTS.T116850766A116851089.enPacoureau, N. et al. Half a century of global decline in oceanic sharks and rays. Nature 589, 567–571 (2021).Article 
ADS 
CAS 

Google Scholar 
Cardeñosa, D. et al. Small fins, large trade: a snapshot of the species composition of low-value shark fins in the Hong Kong markets. Anim. Conserv. 23, 203–211 (2019).Article 

Google Scholar 
Haque, A. B. & Spaet, J. L. Y. Trade in threatened elasmobranchs in the Bay of Bengal, Bangladesh. Fish. Res. 243, 106059 (2021).Article 

Google Scholar 
Alcala, A. C. & Russ, G. R. A direct test of the effects of protective management on abundance and yield of tropical marine resources. ICES J. Mar. Sci. 47, 40–47 (1990).Article 

Google Scholar 
Serrano, A. et al. Effects of anti-trawling artificial reefs on ecological indicators of inner shelf fish and invertebrate communities in the Cantabrian Sea (southern Bay of Biscay). J. Mar. Biol. Assoc. U. Kingd. 91, 623–633 (2011).Article 

Google Scholar 
Cortés, E. Perspectives on the intrinsic rate of population growth. Methods Ecol. Evolution. 7, 1136–1145 (2016).Article 

Google Scholar 
McClenachan, L., Cooper, A. B. & Dulvy, N. K. Rethinking trade-driven extinction risk in marine and terrestrial megafauna. Curr. Biol. 26, 1–7 (2016).Article 

Google Scholar 
Tamburello, N., Cote, I. M. & Dulvy, N. K. Energy and the scaling of animal space use. Am. Naturalist 186, 196–211 (2015).Article 

Google Scholar 
Dulvy, N. K. et al. Challenges and priorities in shark and ray conservation. Curr. Biol. 27, R565–R572 (2017).Article 
CAS 

Google Scholar 
Davidson, L. N. K. & Dulvy, N. K. Global marine protected areas to prevent extinctions. Ecol. Evolution. 1, 1–6 (2017).
Google Scholar 
Pauly, D., Zeller, D. & Palomares, M. L. D. Sea Around Us Concepts, Design and Data, (2021).Simpfendorfer, C. A. & Dulvy, N. K. Bright spots of sustainable shark fishing. Curr. Biol. 27, R83–R102 (2017).Article 

Google Scholar 
Booth, H., Squires, D. & Milner-Gulland, E. J. The mitigation hierarchy for sharks: a risk-based framework for reconciling trade-offs between shark conservation and fisheries objectives. Fish. Fish. 21, 269–289 (2019).Article 

Google Scholar 
Grorud-Colvert, K. et al. The MPA Guide: A framework to achieve global goals for the ocean. Science 373, eabf0861 (2021).Article 
CAS 

Google Scholar 
Enright, S. R., Meneses-Orellana, R. & Keith, I. The Eastern Tropical Pacific Marine Corridor (CMAR): The emergence of a voluntary regional cooperation mechanism for the conservation and sustainable use of marine biodiversity within a fragmented regional ocean governance landscape. Front. Mar. Sci. 8, 674825 (2021).Article 

Google Scholar 
Chin, A., Kyne, P. M., Walker, T. I. & McAuley, R. B. An integrated risk assessment for climate change: analysing the vulnerability of sharks and rays on Australia’s Great Barrier Reef. Glob. Change Biol. 16, 1936–1953 (2010).Article 
ADS 

Google Scholar 
Dwyer, R. G. et al. Individual and population benefits of marine reserves for reef sharks. Curr. Biol. 30, 480–489 (2020).Article 
CAS 

Google Scholar 
Speed, C. W., Cappo, M. & Meekan, M. G. Evidence for rapid recovery of shark populations within a coral reef marine protected area. Biol. Conserv. 220, 308–319 (2018).Article 

Google Scholar 
Mizrahi, M. I., Diedrich, A., Weeks, R. & Pressey, R. L. A systematic review of the socioeconomic factors that influence how marine protected areas impact on ecosystems and livelihoods. Soc. Nat. Resour. 32, 4–20 (2019).Article 

Google Scholar 
IPBES. Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. 1148 (Bonn, Germany, 2019).Butchart, S. H. M. et al. Global biodiversity: indicators of recent declines. Science 328, 1164–1168 (2010).Article 
ADS 
CAS 

Google Scholar 
Hanh, T. T. H. & Boonstra, W. J. What prevents small-scale fishing and aquaculture households from engaging in alternative livelihoods? A case study in the Tam Giang lagoon, Viet Nam. Ocean Coast. Manag. 182, 104943 (2019).Article 

Google Scholar 
Ahmed, N., Troell, M., Allison, E. H. & Muir, J. F. Prawn postlarvae fishing in coastal Bangladesh: challenges for sustainable livelihoods. Mar. Policy. 34, 218–227 (2010).Article 

Google Scholar 
Prasetyo, A. P. et al. Shark and ray trade in and out of Indonesia: addressing knowledge gaps on the path to sustainability. Mar. Policy. 133, 104714 (2021).Article 

Google Scholar 
McClanahan, T., Polunin, N. & Done, T. Ecological states and the resilience of coral reefs. Conserv. Ecol. 6, 18 (2002).
Google Scholar 
Bellwood, D. R., Hughes, T. P. & Hoey, A. S. Sleeping functional group drives coral-reef recovery. Curr. Biol. 16, 2434–2439 (2006).Article 
CAS 

Google Scholar 
Cinner, J. E. et al. Vulnerability of coastal communities to key impacts of climate change on coral reef fisheries. Glob. Environ. Change. 22, 12–20 (2012).Article 
ADS 

Google Scholar 
Víe, J.-C., Hilton-Taylor, C. & Stuart, S. N. Wildlife in a Changing World – An analysis of the 2008 IUCN Red List of Threatened Species. 180 (Gland, Switzerland, 2009).Mace, G. M. et al. Quantification of extinction risk: IUCN’s system for classifying threatened species. Conserv. Biol. 22, 1424–1442 (2008).Article 

Google Scholar 
Sherley, R. B. et al. Estimating IUCN Red List population reduction: JARA – A decision-support tool applied to pelagic sharks. Conserv. Lett. 13, e12688 (2019).
Google Scholar 
IUCN Red List. Threats Classification Scheme (Version 3.2), (2021).Salafsky, N. et al. A standard lexicon for biodiversity conservation: unified classifications of threats and actions. Conserv. Biol. 22, 897–911 (2008).Article 

Google Scholar 
Moore, A. Chiloscyllium arabicum. The IUCN Red List of Threatened Species 2017, e.T161426A109902537 (2017). 10.2305/IUCN.UK.2017-2.RLTS.T161426A109902537.enSadovy de Mitcheson, Y. J. et al. Valuable but vulnerable: Over-fishing and under-management continue to threaten groupers so what now? Mar. Policy. 116, 103909 (2020).Article 

Google Scholar 
R: A language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria, 2021).Regression Models for Ordinal Data v. 2019.12.10 (CRAN, 2019).Econometric Tools for Performance and Risk Analysis v. 2.0.4 (2020).Dormann, C. F. et al. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36, 27–46 (2013).Article 

Google Scholar 
Akinwande, M. O., Dikko, H. G. & Samson, A. Variance inflation factor: As a condition for the inclusion of suppressor variable(s) in regression analysis. Open J. Stat. 5, 754–767 (2015).Article 

Google Scholar 
Burnham, K. P. & Anderson, D. R. Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods Res. 33, 261–304 (2004).Article 
MathSciNet 

Google Scholar 
Plots Coefficients from Fitted Models v. 1.2.8 (2022).Fisheries and Aquaculture Software. FishStatJ – Software for Fishery and Aquaculture Statistical Time Series., (2020).Reynolds, R. W., Rayner, N. A., Smith, T. M., Stokes, D. C. & Wang, W. An improved in situ and satellite SST analysis for climate. J. Clim. 15, 1609–1625 (2002).Article 
ADS 

Google Scholar 
NASA Ocean Biology (OB.DAAC). Mean annual sea surface chlorophyll-a concentration for the period 2009-2013 (composite dataset created by UNEP-WCMC). Data obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua Ocean Colour website (NASA OB.DAAC, Greenbelt, MD, USA), (2014).General Bathymetric Chart of the Oceans. GEBCO_2014 Grid. version 20150318. www.gebco.net (2015).XGBoost: A Scalable Tree Boosting System v. 1.4.1.1 (In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (pp. 785–794). New York, NY, USA: ACM, 2016).Elith, J., Leathwick, J. R. & Hastie, T. A working guide to boosted regression trees. J. Anim. Ecol. 77, 802–813 (2008).Article 
CAS 

Google Scholar 
ArcGIS Pro 2.7.0 (Environmental Systems Research Institute) (2020).Ferreira, L. C. & Simpfendorer, C. Galeocerdo cuvier. The IUCN Red List of Threatened Species 2019, e.T39378A2913541 (2019).Beta Regression v. 3.1-4 (2021).Butchart, S. H. et al. Improvements to the Red List Index. PLoS ONE. 2, e140 (2007).Article 
ADS 

Google Scholar 
Sherman, C. S. et al. Half a century of rising extinction risk of coral reef sharks and rays, sammsherman27/CoralReefSharkRayIUCN: Data and Code Used in Sherman et al. Half a century of rising extinction risk of coral reef sharks and rays v1.0.0. https://doi.org/10.5281/zenodo.7267904 (2022). More

With human movement and globalization, invasive container breeding vectors responsible for dengue, Zika, chikungunya and now malaria, with An. stephensi, are being introduced and establishing populations in new locations. They are bringing with them the threat of increasing or novel cases of vector-borne diseases to new locations where health systems may not be prepared.Anopheles stephensi was first detected on the African continent in Djibouti in 2012 and has since been confirmed in Ethiopia, Somalia, and Sudan. Unlike most malaria vectors, An. stephensi is often found in artificial containers and in urban settings. This unique ecology combined with its initial detection in seaports in Djibouti, Somalia, and Sudan has led scientists to believe that the movement of this vector is likely facilitated through maritime trade.By modeling inter- and intra-continental maritime connectivity in Africa we identified countries with higher likelihood of An. stephensi introduction if facilitated through maritime movement and ranked them based on this data. Anopheles stephensi was not detected in Africa (Djibouti) until 2012. To determine whether historical maritime data would have identified the first sites of introduction, 2011 maritime data were analyzed to determine whether the sites with confirmed An. stephensi would rank highly in connectivity to An. stephensi endemic countries. Using 2011 data on maritime connectivity alone, Djibouti and Sudan were identified as the top two countries at risk of An. stephensi introduction if it is facilitated by marine cargo shipments. In 2021, these are two of the three African coastal nations where An. stephensi is confirmed to be established.When 2011 maritime data were combined with the HSI for An. stephensi establishment, the top five countries remain the same as with maritime data alone: Sudan, Djibouti, Egypt, Kenya and Tanzania, in that order. The maritime data show likelihood of introduction and HSI shows likelihood of establishment. When combined, the analyses show a likelihood of being able to establish and survive once introduced. Interestingly, the results of the combined analyses align with the detection data being reported in the Horn of Africa. The 2011 maritime data reinforces the validity of the model as it points to Sudan and Djibouti, where An. stephensi established in the following years. Similarly, the HSI data for Ethiopia has aligned closely with detections of the species to date15. Interestingly, around this time of initial detection in Djibouti, Djibouti City port underwent development and organizational change. The government of Djibouti took back administrative control of the port as early as 201230.Following this method, maritime trade data from 2020 could point to countries at risk of An. stephensi introduction from endemic countries as well as from the coastal African countries with newly introduced populations. Here we provide a prioritization list and heat map of countries for the early detection, rapid response, and targeted surveillance of An. stephensi in Africa based on this data and the HSI (Fig. 4). Further invasion of An. stephensi on the African continent has the potential to reverse progress made on malaria control in the last century. Anopheles stephensi thrives in urban settings and in containers, in contrast to the rural settings and natural habitats where most Anopheles spp. are found20. The situation in Djibouti may be a harbinger for what is to come if immediate surveillance and control strategies are not initiated18.Figure 4Prioritization Heat Map of African Countries. These 2020 heat maps rank African countries using (A) the Likelihood of An. stephensi through Maritime Trade Index (LASIMTI) data alone and (B) LASIMTI and HSI combined, based on maritime connectivity to countries where An. stephensi is endemic. Higher ranking countries which are at greater risk of An. stephensi introduction are darker in red color than those that are lower ranking (lighter red). Countries which are shaded grey are inland countries that do not have a coast and therefore no data on maritime movement into ports. Countries which are grey and checkered have established or endemic An. stephensi populations and are considered source locations for potential An. stephensi introduction in this analysis. Map was generated using MapChart (mapchart.net).Full size imageMaritime data from 2020, with Djibouti and Sudan considered as potential source populations for intracontinental introduction of An. stephensi, indicate the top five countries at risk for maritime introduction are Egypt, Kenya, Mauritius, Tanzania, and Morocco, suggesting that targeted larval surveillance in these countries near seaports may provide a better understanding of whether there are maritime introductions. When the data from 2020 data is combined with HSI for An. stephensi, the top five countries are instead Egypt, Kenya, Tanzania, Morocco, and Libya. Interestingly, historical reports of An. stephensi in Egypt exist; however, following further identification these specimens were determined to be An. ainshamsi31. With several suitable habitats both along the coast and inland of Egypt, revisiting surveillance efforts there would provide insight into how countries that are highly connected to An. stephensi locations through maritime traffic may experience introductions.Further field validation of this prioritization list is necessary, because it is possible that An. stephensi is being introduced through other transportation routes, such as dry ports or airports32, or may even be dispersed through wind facilitation33. However, countries highlighted here with high levels of connectivity to known An. stephensi locations should be considered seriously at risk and surveillance urgently established to determine whether An. stephensi introduction has already occurred or to enable early detection. Primary vector surveillance for both Ae. aegypti and An. stephensi are through larval surveys, and the two mosquitoes are commonly detected in the same breeding habitats. It could therefore be beneficial to coordinate with existing Aedes surveillance efforts to be able to simultaneously gather data on medically relevant Aedes vectors while seeking to determine whether An. stephensi is present. Similarly, in locations with known An. stephensi and not well established Aedes programs, coordinating surveillance efforts provides an opportunity to conduct malaria and arboviral surveillance by container breeding mosquitoes simultaneously.Efforts to map pinch points or key points of introduction based on the movement of goods and populations could provide high specificity for targeted surveillance and control efforts. For example, participatory mapping or population mobility data collection methods, such as those used to determine routes of human movement for malaria elimination, may simultaneously provide information on where targeted An. stephensi surveillance efforts should focus. Several methods have been proposed in the literature for modeling human movement and one in particular, PopCAB, which is often used for communicable diseases, combined quantitative and qualitative data with geospatial information to identify points of control34.Data on invasive mosquito species has shown that introduction events are rarely a one-time occurrence. Population genetics data on Aedes species indicate that reintroductions are very common and can facilitate the movement of genes between geographically distinct populations, raising the potential for introduction of insecticide resistance, thermotolerance, and other phenotypic and even behavioral traits which may be facilitated by gene flow and introgression35. Djibouti, Sudan, Somalia, and Ethiopia, countries with established invasive populations of An. stephensi, should continue to monitor invasive populations and points of introduction to control and limit further expansion and adaptation of An. stephensi. Work by Carter et al. has shown that An. stephensi populations in Ethiopia in the north and central regions can be differentiated genetically, potentially indicating that these populations are a result of more than one introduction into Ethiopia from South Asia, further emphasizing the potential role of anthropogenic movement on the introduction of the species17.One major limitation of this work is that Somalia is the third coastal nation where An. stephensi has been confirmed; however, marine traffic data were not available for Somalia so it could not be included in this analysis. The potential impact of Somalia on maritime trade is unknown and it should not be excluded as a potential source population. Additionally, this model does not account for the possibility of other countries with An. stephensi populations that have not been detected yet. As new data on An. stephensi expansion becomes available, more countries will be at higher risk. Other countries with An. stephensi populations, such as Iran, Myanmar, and Iraq, constitute lower relative percentages of trade with these countries so were not included in the analysis. However, genetic similarities were noted from An. stephensi in Pakistan, so this nation was included10.Due to the nature of maritime traffic, inland countries were also not included in this prioritization ranking. Countries which are inland but share borders with high-risk countries according to the LASTIMI index should also be considered with high priority. For example, the ranking from 2011 highlights Sudan and Djibouti, both which border Ethiopia, and efforts to examine key land transportation routes between bordering nations where humans and goods travel may provide additional insight into the expansion routes of this invasive species.In Ethiopia, An. stephensi was detected in 2016. It has largely been detected along major transportation routes although further data is needed to understand the association between movement and An. stephensi introductions and expansion since most sampling sites have also been located along transport routes. Importantly, Ethiopia relies heavily on the ports of Djibouti and Somalia for maritime imports and exports. Surveillance efforts have revealed that the species is also frequently associated with livestock shelters and An. stephensi are frequently found with livestock bloodmeals15. Interestingly, the original detection of An. stephensi was found in a livestock quarantine station in the port of Djibouti. Additionally, livestock constitutes one of the largest exports of maritime trade from this region. For countries with high maritime connectivity to An. stephensi locations, surveillance efforts near seaports, in particular those with livestock trade, may be targeted locations for countries without confirmed An. stephensi to begin larval surveillance.As Ae. aegypti and Culex coronator were detected in tires or Ae. albopictus through tire and bamboo (Dracaena sanderiana) trade, An. stephensi could be carried through maritime trade of a specific good36,37,38. Future examination of the movement of specific goods would be beneficial in interpreting potential An. stephensi invasion pathways. Additionally, the various types of vessels used to transport certain cargo such as container, bulk, and livestock ships could affect An. stephensi survivability during transit. Sugar and grain are often shipped in bulk or break bulk vessels which store cargo in large unpackaged containers. Container ships transport products stored in containers sized for land transportation via trucks and carry goods such as tires. Livestock vessels are often multilevel, open-air ships which require more hands working on deck and water management39.Using LSBCI index data from 2020, we developed a network to highlight how coastal African nations are connected through maritime trade (Fig. 4). The role of this network analysis is two-fold, (1) it demonstrates an understanding of intracontinental maritime connectivity; and (2) it highlights the top three countries connected via maritime trade through an interactive html model (Supplemental File). For example, if An. stephensi is detected and established in a specific coastal African nation such as Djibouti, selecting the Djibouti node reveals the top three locations at risk of introduction from that source country (Djibouti links to Sudan, Egypt and Kenya). This can be used as an actionable prioritization list for surveillance if An. stephensi is detected in any given country and highlights major maritime hubs in Africa which could be targeted for surveillance and control. For example, since the development of this model, An. stephensi has been detected in Nigeria. Through the use of this interactive model, Ghana, Cote d’Ivoire, and Benin have been identified as countries most connected to Nigeria through maritime trade and therefore surveillance prioritization activities could consider these locations.The network analysis reveals the significance of South African trade to the rest of the continent. Due to the distance, South Africa did not appear to be high in risk of An. stephensi introduction. However, this analysis does reveal that if An. stephensi were to enter nearby countries, it could very easily be introduced because of its high centrality. Western African countries such as Ghana, Togo, and Morocco are also heavily connected to other parts of Africa. Interestingly, Mauritius appears to be highly significant to this network of African maritime trade. Based on 2020 maritime data, Mauritius is ranked as the country with the third greatest likelihood of introduction of An. stephensi and has the second highest centrality rank value of 0.159. Considering these factors, Mauritius could serve as an important port of call connecting larger ports throughout Africa or other continents. With long standing regular larval surveillance efforts across the island for Aedes spp., this island nation is well suited to look for Anopheles larvae as part of Aedes surveillance efforts for early detection and rapid response to prevent the establishment of An. stephensi. If An. stephensi were to become established in countries with high centrality ranks, further expansion on the continent could be accelerated drastically. These ports could serve as important watchpoints and indicators of An. stephensi’s incursion into Africa. Anopheles stephensi is often found in shared habitats with Aedes spp. and a great opportunity exists to leverage Aedes arboviral surveillance efforts to initiate the search for An. stephensi, especially in countries that have high potential of introduction through maritime trade. More

In 2008, I was investigating the methane bubbling up on the beaches and in shallow waters of Mocha Island, off the coast of central Chile. I became intrigued by how microorganisms could thrive in methane-rich areas and changed my research focus from marine biology to extreme environments. I wanted to understand how methane acts as a source of energy and carbon for microbes.Since then, I have explored a number of bizarre environments. In 2010, I went in a submarine down to 200 metres in the Black Sea, one of the world’s largest anoxic water bodies. There, I found mats of filamentous bacteria that survive on sulfur compounds.In 2017, I studied the microbes in Canada’s tailing ponds, artificial lakes of water, sand and clay waste that are left behind after petroleum extraction. And I sampled the microorganisms living in 100 °C Antarctic hot springs in 2022.I came home to Chile in 2018 and began collaborating with an international team researching the geomicrobiology of thermal features, including hot springs, geysers and volcanoes. After travelling with the group to Argentina’s active volcanic region, I got funding to explore the microbial communities that exist beneath hydrothermal vents in southern Chile, where the oceanic crust is subducting beneath the continental plate.In this image, I am in the Atacama Desert in South America, the driest non-polar desert on the planet. I am measuring 80–100 °C steam released from a fumarole containing yellow sulfur, which crystallizes at its opening as the vapour cools. I also sampled sub-surface microbes that are flushed out with the fluids. We’ll sequence their DNA to assess the microbial communities and their biological interactions.My goal is to learn more about subsurface microbes in extreme environments. I want to understand how microbial forces shaped the planet and how these communities might shift in the future with climate change. More

Stefansson, S. O., Björnsson, B. T., Ebbesson, L. O. E. & McCormick, S. D. Smoltification. In Fish Larval Physiology (eds Finn, R. N. & Kapoor, B. G.) 639–681 (CRC Press, 2020).Chapter 

Google Scholar 
Kaleka, A. S., Kaur, N. & Bali, G. K. Larval development and molting. In Edible Insects (ed. Mikkola, H.) 17 (IntechOpen, 2019).
Google Scholar 
Butler, L. K. & Rohwer, V. G. Feathers and molt. in Ornithology: Foundation, Analysis, and Application (eds Morrison, M. L. et al.) 242–270 (JHU Press, 2018).
Google Scholar 
Swaddle, J. P., Witter, M. S., Cuthill, I. C., Budden, A. & McCowen, P. Plumage condition affects flight performance in common starlings: Implications for developmental homeostasis, abrasion and moult. J. Avian Biol. 27, 103–111 (1996).Article 

Google Scholar 
Norris, D. R., Marra, P. P., Montgomerie, R., Kyser, T. K. & Ratcliffe, L. M. Reproductive effort, molting latitude, and feather color in a migratory songbird. Science 306, 2249–2250 (2004).Article 
ADS 
CAS 

Google Scholar 
Delhey, K., Peters, A. & Kempenaers, B. Cosmetic coloration in birds: Occurrence, function, and evolution. Am. Nat. 169, S145–S158 (2007).Article 

Google Scholar 
Tomotani, B. M. & Muijres, F. T. A songbird compensates for wing molt during escape flights by reducing the molt gap and increasing angle of attack. J. Exp. Biol. 222, 195396 (2019).Article 

Google Scholar 
Galván, I., Negro, J. J., Rodriguez, A. & Carrascal, L. M. On showy dwarfs and sober giants: Body size as a constraint for the evolution of bird plumage colouration. Acta Ornithol. 48, 65–80 (2013).Article 

Google Scholar 
Speakman, J. R. & Król, E. Maximal heat dissipation capacity and hyperthermia risk: Neglected key factors in the ecology of endotherms. J. Anim. Ecol. 79, 726–746 (2010).
Google Scholar 
Wolf, B. O. & Walsberg, G. E. The role of the plumage in heat transfer processes of birds. Am. Zool. 40, 575–584 (2000).
Google Scholar 
Berthold, P. & Querner, U. Genetic basis of moult, wing length, and body weight in a migratory bird species, Sylvia atricapilla. Experientia 38, 801–802 (1982).Article 

Google Scholar 
Gwinner, E., Neusser, V., Engl, D., Schmidl, D. & Bals, L. Haltung, Zucht und Eiaufzucht afrikanischer und europäischer Schwarzkehlchen Saxicola torquata. Gefied. Welt 111, 118–120 (1987).
Google Scholar 
Berthold, P. & Querner, U. Microevolutionary aspects of bird migration based on experimental results. Isr. J. Ecol. Evol. 41, 377–385 (1995).
Google Scholar 
Helm, B. & Gwinner, E. Timing of postjuvenal molt in African (Saxicola torquata axillaris) and European (Saxicola torquata rubicola) stonechats: Effects of genetic and environmental factors. Auk 116, 589–603 (1999).Article 

Google Scholar 
Helm, B. & Gwinner, E. Timing of molt as a buffer in the avian annual cycle. Acta Zool. Sin. 52, 703–706 (2006).
Google Scholar 
Rohwer, S., Ricklefs, R. E., Rohwer, V. G. & Copple, M. M. Allometry of the duration of flight feather molt in birds. PLoS Biol. 7, e1000132 (2009).Article 

Google Scholar 
Jenni, L. & Winkler, R. The Biology of Moult in Birds (Bloomsbury Publishing, 2020).
Google Scholar 
Tonra, C. M. & Reudink, M. W. Expanding the traditional definition of molt-migration. Auk Ornithol. Adv. 135, 1123–1132 (2018).
Google Scholar 
Rohwer, S., Butler, L. K., Froehlich, D. R., Greenberg, R. & Marra, P. P. Ecology and demography of east–west differences in molt scheduling of Neotropical migrant passerines. Birds Two Worlds Ecol. Evol. Migr. (R. Greenb. PP Marra, Eds.). Johns Hopkins Univ. Press. Balt. Maryl., 87–105 (2005).Bensch, S., Åkesson, S. & Irwin, D. E. The use of AFLP to find an informative SNP: Genetic differences across a migratory divide in willow warblers. Mol. Ecol. 11, 2359–2366 (2002).Article 
CAS 

Google Scholar 
Ruegg, K. Genetic, morphological, and ecological characterization of a hybrid zone that spans a migratory divide. Evol. Int. J. Org. Evol. 62, 452–466 (2008).Article 

Google Scholar 
Delmore, K. E., Fox, J. W. & Irwin, D. E. Dramatic intraspecific differences in migratory routes, stopover sites and wintering areas, revealed using light-level geolocators. Proc. R. Soc. B Biol. Sci. 279, 4582–4589 (2012).Article 

Google Scholar 
Delmore, K. E. et al. Individual variability and versatility in an eco-evolutionary model of avian migration. Proc. R. Soc. B 287, 20201339 (2020).Article 

Google Scholar 
Procházka, P. et al. Across a migratory divide: divergent migration directions and non-breeding grounds of Eurasian reed warblers revealed by geolocators and stable isotopes. J. Avian Biol. 49, 012516 (2018).Article 

Google Scholar 
Bensch, S., Grahn, M., Müller, N., Gay, L. & Åkesson, S. Genetic, morphological, and feather isotope variation of migratory willow warblers show gradual divergence in a ring. Mol. Ecol. 18, 3087–3096 (2009).Article 

Google Scholar 
Rohwer, S. & Irwin, D. E. Molt, orientation, and avian speciation. Auk 128, 419–425 (2011).Article 

Google Scholar 
Pageau, C., Sonnleitner, J., Tonra, C. M., Shaikh, M. & Reudink, M. W. Evolution of winter molting strategies in European and North American migratory passerines. Ecol. Evol. 11, 13247–13258 (2021).Article 

Google Scholar 
Butler, L. K., Rohwer, S. & Rogers, M. Prebasic molt and molt-related movements in Ash-throated Flycatchers. Condor 108, 647–660 (2006).Article 

Google Scholar 
Barry, J. H., Butler, L. K., Rohwer, S. & Rohwer, V. G. Documenting molt-migration in Western Kingbird (Tyrannus verticalis) using two measures of collecting effort. Auk 126, 260–267 (2009).Article 

Google Scholar 
Hobson, K. A. & Wassenaar, L. I. Linking breeding and wintering grounds of neotropical migrant songbirds using stable hydrogen isotopic analysis of feathers. Oecologia 109, 142–148 (1996).Article 
ADS 
CAS 

Google Scholar 
Hobson, K. A. & Wassenaar, L. I. Tracking Animal Migration with Stable Isotopes (Academic Press, 2018).
Google Scholar 
Rubenstein, D. R. & Hobson, K. A. From birds to butterflies: Animal movement patterns and stable isotopes. Trends Ecol. Evol. 19, 256–263 (2004).Article 

Google Scholar 
Bearhop, S. et al. Assortative mating as a mechanism for rapid evolution of a migratory divide. Science 310, 502–504 (2005).Article 
ADS 
CAS 

Google Scholar 
Eppig, J. T. et al. The mouse genome database (MGD): Comprehensive resource for genetics and genomics of the laboratory mouse. Nucleic Acids Res. 40, D881–D886 (2012).Article 
CAS 

Google Scholar 
Contina, A., Bridge, E. S. & Kelly, J. F. Exploring novel candidate genes from the mouse genome informatics database: Potential implications for avian migration research. Integr. Zool. 11, 240 (2016).Article 

Google Scholar 
Yang, J. et al. Common SNPs explain a large proportion of the heritability for human height. Nat. Genet. 42, 565–569 (2010).Article 
CAS 

Google Scholar 
Thompson, C. W. Is the Painted Bunting actually two species? Problems determining species limits between allopatric populations. Condor 93, 987–1000 (1991).Article 

Google Scholar 
Contina, A., Bridge, E. S., Seavy, N. E., Duckles, J. M. & Kelly, J. F. Using geologgers to investigate bimodal isotope patterns in Painted Buntings (Passerina ciris). Auk 130, 265 (2013).Article 

Google Scholar 
Besozzi, E., Chew, B., Allen, D. C. & Contina, A. Stable isotope analysis of an aberrant Painted Bunting (Passerina ciris) feather suggests post-molt movements. Wilson J. Ornithol. 133, 151 (2021).Article 

Google Scholar 
Sharp, A. et al. Spatial and Temporal Scale-Dependence of the Strength of Migratory Connectivity in a North American Passerine. https://assets.researchsquare.com/files/rs-1483049/v1/72236b63-952d-4870-89e7-461056b8625b.pdf?c=1648893558 (2022).Pyle, P. et al. Temporal, spatial, and annual variation in the occurrence of molt-migrant passerines in the Mexican monsoon region. Condor 111, 583–590 (2009).Article 

Google Scholar 
Bridge, E. S., Fudickar, A. M., Kelly, J. F., Contina, A. & Rohwer, S. Causes of bimodal stable isotope signatures in the feathers of a molt-migrant songbird. Can. J. Zool. 89, 951 (2011).Article 
CAS 

Google Scholar 
Seutin, G., White, B. N. & Boag, P. T. Preservation of avian blood and tissue samples for DNA analyses. Can. J. Zool. 69, 82–90 (1991).Article 
CAS 

Google Scholar 
Ali, O. A. et al. RAD capture rapture: Flexible and efficient sequence-based genotyping. Genetics 202, 389–400 (2016).Article 
CAS 

Google Scholar 
Contina, A. et al. Characterization of SNP markers for the Painted Bunting (Passerina ciris) and their relevance in population differentiation and genome evolution studies. Conserv. Genet. Resour. 11, 5–10 (2019).Article 
ADS 

Google Scholar 
Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A. & Cresko, W. A. Stacks: An analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140 (2013).Article 

Google Scholar 
Parker, P., Li, B., Li, H. & Wang, J. The genome of Darwin’s Finch (Geospiza fortis). Gigascience 10, 100040 (2012).
Google Scholar 
Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012).Article 
CAS 

Google Scholar 
Van der Auwera, G. A. et al. From FastQ data to high-confidence variant calls: The genome analysis toolkit best practices pipeline. Curr. Protoc. Bioinform. 43, 1–33 (2013).
Google Scholar 
McKenna, A. et al. The genome analysis toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010).Article 
CAS 

Google Scholar 
Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156–2158 (2011).Article 
CAS 

Google Scholar 
Anderson, E. genoscapeRtools: Tools for Building Migratory Bird Genoscapes (2019).Purcell, S. et al. PLINK: A tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007).Article 
CAS 

Google Scholar 
Alexander, D. H., Novembre, J. & Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664 (2009).Article 
CAS 

Google Scholar 
Alexander, D. H. & Lange, K. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinform. 12, 246 (2011).Article 

Google Scholar 
Francis, R. M. pophelper: An R package and web app to analyse and visualize population structure. Mol. Ecol. Resour. 17, 27–32 (2017).Article 
CAS 

Google Scholar 
Chew, B., Kelly, J. & Contina, A. Stable isotopes in avian research: a step by step protocol to feather sample preparation for stable isotope analysis of carbon (δ13C), nitrogen (δ15N), and hydrogen (δ2H). Version 1.1. https://doi.org/10.17504/protocols.io.z2uf8ew (2019).Wassenaar, L. I. & Hobson, K. A. Comparative equilibration and online technique for determination of non-exchangeable hydrogen of keratins for use in animal migration studies. Isotopes Environ. Health Stud. 39(3), 211–217 (2003).Article 
CAS 

Google Scholar 
Bowen, G. J., Wassenaar, L. I. & Hobson, K. A. Global application of stable hydrogen and oxygen isotopes to wildlife forensics. Oecologia 143, 337–348 (2005).Article 
ADS 

Google Scholar 
R Core Team: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2021).Wassenaar, L. I. & Hobson, K. A. Stable-hydrogen isotope heterogeneity in keratinous materials: Mass spectrometry and migratory wildlife tissue subsampling strategies. Rapid Commun. Mass Spectrom. 20, 2505–2510 (2006).Article 
ADS 
CAS 

Google Scholar 
Zhou, X., Carbonetto, P. & Stephens, M. Polygenic modeling with Bayesian sparse linear mixed models. PLoS Genet. 9, e1003264 (2013).Article 
CAS 

Google Scholar 
Guan, Y. & Stephens, M. Bayesian variable selection regression for genome-wide association studies and other large-scale problems. Ann. Appl. Stat. 5, 455 (2011).Article 
MathSciNet 
MATH 

Google Scholar 
Marchini, J., Cardon, L. R., Phillips, M. S. & Donnelly, P. The effects of human population structure on large genetic association studies. Nat. Genet. 36, 512–517 (2004).Article 
CAS 

Google Scholar 
Browning, S. R. & Browning, B. L. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. Am. J. Hum. Genet. 81, 1084–1097 (2007).Article 
CAS 

Google Scholar 
Chaves, J. A. et al. Genomic variation at the tips of the adaptive radiation of Darwin’s finches. Mol. Ecol. 25, 5282–5295 (2016).Article 
CAS 

Google Scholar 
Zhang, Y.-W. et al. mrMLM v4.0.2: An R platform for multi-locus genome-wide association studies. Genom. Proteom. Bioinform. 18, 481–487 (2020).Article 

Google Scholar 
Grabherr, M. G. et al. Genome-wide synteny through highly sensitive sequence alignment: Satsuma. Bioinformatics 26, 1145–1151 (2010).Article 
CAS 

Google Scholar 
Quinlan, A. R. & Hall, I. M. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842 (2010).Article 
CAS 

Google Scholar 
Ellis, N., Smith, S. J. & Pitcher, C. R. Gradient forests: Calculating importance gradients on physical predictors. Ecology 93, 156–168 (2012).Article 

Google Scholar 
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).Article 

Google Scholar 
Anderson, E. C. snps2assays: Prepare SNP Assay Orders from ddRAD or RAD Loci (2015).Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760 (2009).Article 
CAS 

Google Scholar 
Patterson, N., Price, A. L. & Reich, D. Population structure and eigenanalysis. PLoS Genet. 2, e190 (2006).Article 

Google Scholar 
Price, A. L. et al. Principal components analysis corrects for stratification in genome-wide association studies. Nat. Genet. 38, 904–909 (2006).Article 
CAS 

Google Scholar 
Ruegg, K. et al. Ecological genomics predicts climate vulnerability in an endangered southwestern songbird. Ecol. Lett. 21, 1085–1096 (2018).Article 

Google Scholar 
Bay, R. A. et al. Genomic signals of selection predict climate-driven population declines in a migratory bird. Science 359, 83–86 (2018).Article 
ADS 
CAS 

Google Scholar 
Hedenström, A. Adaptations to migration in birds: Behavioural strategies, morphology and scaling effects. Philos. Trans. R. Soc. B Biol. Sci. 363, 287–299 (2008).Article 

Google Scholar 
Buehler, D. M. & Piersma, T. Travelling on a budget: Predictions and ecological evidence for bottlenecks in the annual cycle of long-distance migrants. Philos. Trans. R. Soc. B Biol. Sci. 363, 247–266 (2008).Article 

Google Scholar 
Schieltz, P. C. & Murphy, M. E. The contribution of insulation changes to the energy cost of avian molt. Can. J. Zool. 75, 396–400 (1997).Article 

Google Scholar 
Carling, M. D. & Thomassen, H. A. The role of environmental heterogeneity in maintaining reproductive isolation between hybridizing Passerina (Aves: Cardinalidae) buntings. Int. J. Ecol. 2012, 1–11 (2012).Article 

Google Scholar 
Irwin, D. E. Incipient ring speciation revealed by a migratory divide. Mol. Ecol. 18, 2923–2925 (2009).Article 

Google Scholar 
Thomas, D. W., Blondel, J., Perret, P., Lambrechts, M. M. & Speakman, J. R. Energetic and fitness costs of mismatching resource supply and demand in seasonally breeding birds. Science 291, 2598–2600 (2001).Article 
ADS 
CAS 

Google Scholar 
Rohwer, V. G., Rohwer, S. & Ortiz-Ramirez, M. F. Molt biology of resident and migrant birds of the monsoon region of west Mexico. Ornitol. Neotrop. 20, 565–584 (2009).
Google Scholar 
Bensch, S., Andersson, T. & Åkesson, S. Morphological and molecular variation across a migratory divide in willow warblers, Phylloscopus trochilus. Evolution 53, 1925–1935 (1999).Article 

Google Scholar 
Turbek, S. P., Scordato, E. S. C. & Safran, R. J. The role of seasonal migration in population divergence and reproductive isolation. Trends Ecol. Evol. 33, 164–175 (2018).Article 

Google Scholar 
Scordato, E. S. C. et al. Migratory divides coincide with reproductive barriers across replicated avian hybrid zones above the Tibetan Plateau. Ecol. Lett. 23, 231–241 (2020).Article 

Google Scholar 
Battey, C. J. et al. A migratory divide in the Painted Bunting (Passerina ciris). Am. Nat. 191, 259–268 (2018).Article 
CAS 

Google Scholar 
Contina, A. et al. Genetic structure of the Painted Bunting and its implications for conservation of migratory populations. Ibis 161, 372 (2019).Article 

Google Scholar 
Butler, L. K. The grass is always greener: Do monsoon rains matter for molt of the Vermilion Flycatcher (Pyrocephalus rubinus)? Auk 130, 297–307 (2013).Article 

Google Scholar 
Turbek, S. P. et al. A migratory divide spanning two continents is associated with genomic and ecological divergence. Evolution 76, 722 (2022).Article 

Google Scholar 
Dietz, M. W., Daan, S. & Masman, D. Energy requirements for molt in the kestrel Falco tinnunculus. Physiol. Zool. 65, 1217–1235 (1992).Article 

Google Scholar 
Vézina, F., Gustowska, A., Jalvingh, K. M., Chastel, O. & Piersma, T. Hormonal correlates and thermoregulatory consequences of molting on metabolic rate in a northerly wintering shorebird. Physiol. Biochem. Zool. 82, 129–142 (2009).Article 

Google Scholar 
Bazzi, G. et al. Candidate genes have sex-specific effects on timing of spring migration and moult speed in a long-distance migratory bird. Curr. Zool. 63, 479–486 (2017).CAS 

Google Scholar 
Busby, L. et al. Sonic hedgehog specifies flight feather positional information in avian wings. Development 147, 188821 (2020).Article 

Google Scholar 
Eichberger, T. et al. GLI2-specific transcriptional activation of the bone morphogenetic protein/Activin antagonist Follistatin in human epidermal cells. J. Biol. Chem. 283, 12426–12437 (2008).Article 
CAS 

Google Scholar 
Matzuk, M. M. et al. Multiple defects and perinatal death in mice deficient in follistatin. Nature 374, 360–363 (1995).Article 
ADS 
CAS 

Google Scholar 
Patel, K., Makarenkova, H. & Jung, H.-S. The role of long range, local and direct signalling molecules during chick feather bud development involving the BMPs, follistatin and the Eph receptor tyrosine kinase Eph-A4. Mech. Dev. 86, 51–62 (1999).Article 
CAS 

Google Scholar 
Nakamura, M. et al. Control of pelage hair follicle development and cycling by complex interactions between follistatin and activin. FASEB J. 17, 1–22 (2003).Article 
MathSciNet 

Google Scholar 
Pays, L., Charvet, I., Hemming, F. J. & Saxod, R. Close link between cutaneous nerve pattern development and feather morphogenesis demonstrated by experimental production of neo-apteria and ectopic feathers: Implication of chondroitin sulphate proteoglycans and other matrix molecules. Anat. Embryol. 195, 457–466 (1997).Article 
CAS 

Google Scholar 
Pyle, P., Saracco, J. F. & DeSante, D. F. Evidence of widespread movements from breeding to molting grounds by North American landbirds. Auk Ornithol. Adv. 135, 506–520 (2018).
Google Scholar 
De Mita, S. et al. Detecting selection along environmental gradients: Analysis of eight methods and their effectiveness for outbreeding and selfing populations. Mol. Ecol. 22, 1383–1399 (2013).Article 

Google Scholar 
Lotterhos, K. E. & Whitlock, M. C. Evaluation of demographic history and neutral parameterization on the performance of FST outlier tests. Mol. Ecol. 23, 2178–2192 (2014).Article 

Google Scholar 
Frichot, E., Schoville, S. D., de Villemereuil, P., Gaggiotti, O. E. & François, O. Detecting adaptive evolution based on association with ecological gradients: Orientation matters!. Heredity (Edinb.) 115, 22–28 (2015).Article 
CAS 

Google Scholar 
Trivedi, A. K. et al. Temperature alters the hypothalamic transcription of photoperiod responsive genes in induction of seasonal response in migratory redheaded buntings. Mol. Cell. Endocrinol. 493, 110454 (2019).Article 
CAS 

Google Scholar  More

Morphometric measurements are subjective and used to assess the carcass development and quantitatively measure the muscular distribution in the carcass with estimates of its conformation. In the present study there were not significative differences observed for these parameters or for carcass compactness index (CCI), inferring that the use of cactus pear silage as well as intermittent water supply combined or alone did not alter animal growth and/or carcass conformation, maintaining the muscle pattern achieved by the control diet (usual) and demonstrating body and carcass uniformity. Since animals used in this study were homogeneous and had similar age and body performance, as indicated by the carcass morphometric measurements and by the difference between the empty carcass and hot carcass weights, which resulted in the sum of head + limb with an average of 8.2 ± 0.13 kg between treatments, giving an idea that the animals were similar in chronological age, since the allometric growth of the body occurs from the extremities to the interior of the body.The significant difference between treatments with inclusion of cactus pear silage for hot carcass yield (HCY) and cold carcass yield (CCY) may be related to the weight of the full gastrointestinal tract, which showed higher values for animals fed with a higher proportion of Tifton 85 grass hay in the diet (0% CPS). Increasing the NDF content of the diet reduces the passage rate of digesta, and the emptying of the gastrointestinal tract (GT) that cause a distension of the rumen-reticulum and increase the weight of the gastrointestinal tract, resulting in lower HCY and consequently lower CCY. While the diets with inclusion of CPS increase NFC content, such as pectin, which have higher rates of rumen degradability and, higher rates of passage7,8,9.Measurements and evaluations carried out on the carcass, such as the carcass compactness index and loin eye area (LEA), are parameters that quantitatively measure the muscle distribution in the carcass, an edible part of greater financial return, which indicates the conformation of these animals3, while the body condition score (BCS) and the measure C, which are highly correlated, measure the distribution of fat on the carcass, giving an idea of the carcass finish, in which the higher these variables, the greater the proportion of fat that allows for less water loss due to carcass cooling10. These variables in the present study were also not influenced by the levels of cactus pear silage and water restrictions, presenting an overall mean of 0.17 kg/cm, 7.68 cm, 2.42 points and 0.7 mm respectively, and consequently did not influence the losses due to cooling, which presented an average loss of 1.48%.The main cuts of the goat carcass are the neck, leg, shoulder, loin, and rib. Their economic values differ, and their proportions become an important index to evaluate the carcass quality9. The cuts of greatest importance and commercial values are the leg and the loin, called noble cuts because they present greater yield and muscle tenderness, being interesting that they present a good proportion in the carcass, for providing greater edible tissue content, mainly muscle.Carcasses with similar weight tend to have equivalent proportions of cuts, as they exhibit isogonic growth. As the cold carcass weight (CCW) and the conformation of the animals were similar, with similar morphometric measurements, they had a direct relationship in the absence of an effect on commercial cuts.The commercial value of the carcass, whether through carcass yield and/or the proportions of the cuts, is also linked to tissue composition, thus the dissection of the leg represents an estimate of measuring the tissue composition of the carcass, in which is sought a greater proportion of muscle, intermediate proportion of fat and less bone in carcasses11. In this way, diets with cactus pear silage and the different levels of intermittent water supply resulted in the constancy in the amount of muscle, fat, and bone in legs of goats. The similarity in muscle proportion is related to the lack of effects on slaughter weight and CCW, as the weight of muscles is highly correlated to carcass weight. The average muscle yield was above 60% in all treatments, confirming that the animals showed good efficiency to the diets and adapted well to the water supply levels. Although the diets with cactus silage had high amounts of metabolizable energy (ME) and no difference in DM intake, the energy input was similar that not influencing carcass weights and carcass compactness index. That is, it did not influence muscle deposition in the carcass, probably due to synchronicity of energy and protein.As for the weight and proportion of bone tissue, it is believed that because this is a tissue with early development in relation to muscle and fat2, diets in the final stages of growth (average of 8 months) would hardly change their participation in the tissue composition, where the relationship of this tissue with the others is usually only increased when there are changes in the proportion of muscle and/or fat.Water restriction, as long as it is moderate and acute, mainly affects the loss of body water and not tissues, which does not cause deleterious effects on animal productivity and growth.The muscle:fat ratio indicates the state of leg fattening, while the muscle:bone ratio estimates the carcass muscularity, both being attributes of quality3. The similarity previously reported in the weight of fat, bone and muscle corroborates that these relationships also do not have differences. The same occurs for the leg muscularity index (LMI), due to the weight of the five muscles used to determine the index and the length of the femur which had been similar between the animals.Nevertheless, when considering fat as a percentage of participation in leg weight, it is possible to observe that the intermittency in water supply in both intervals (24 and 48 h) reduced the proportion of fat in the leg. Although in this research, the water supply levels did not affect the daily intake of dry matter from animals, with average intake of 650.67 g/kg DM, ranging from 599 to 682 g/kg DM between treatments7, during days of water deprivation, fat mobilization for energy availability may occur, possibly offsetting water stress and influencing not only feed intake, on these days of deprivation but also affecting energy metabolism, which results in the mobilization of energy reserves2.When the physicochemical composition of the meat was evaluated, it was observed that the diets and water supply levels probably did not affect the reserves of muscle glycogen during the pre-slaughter management as can be seen through pHinitial and pHfinal. The pHinitial right after slaughter should be close to neutrality, as well as in the live animal, indicating that the animal did not suffer from stress during the pre-slaughter period. The pHfinal, on the other hand, is expected to show a considerable variation, between 5.55 and 6.2 for goat meat; and due be inversely proportional to the concentration of muscle glycogen at the time of slaughter, that is, a more intense expenditure of glycogen stores results in less lactic acid production and higher pHfinal10,12,13. In this research, the pHfinal had an average of 5.74, a pH higher than the isoelectric point of muscle proteins (5.2–5.3). This result is favorable, since it is above the neutral charge and presenting an excessive negative charge that provides the repulsion of filaments, which allows water molecules to bind and improve the organoleptic characteristics of the meat, through succulence and texture of meat13 evaluated by cooking loss, moisture, and shear force, principally. The cooking loss (CL), moisture and shear force (SF) were within the values recommended (20–35% CL, moisture above 70% and SF up to 44.13 Newton (N) for goat meat) to classify the meat as soft and tender14. Statistically, interactions were found between the supply of silage and intermittent water supply, in which goats on a diet without cactus pear silage and without intermittent water supply showed higher values of cooking losses and shear force.Higher concentrations of collagen content and/or greater activities of calpastatin (which inhibit the action of calpains), as well as larger fascicles and greater number of fibers present in each muscle fascicle, as was visually observed in the meat of the animals in this research, can lead to reductions in meat tenderness15. Because goat carcasses are generally small, with low marbling degree and a thin layer of subcutaneous fat, there is rapid heat dissipation at the beginning of the post-mortem period, which can lead to cold shortening, muscle hardening, and less tender meats16.pHfinal of the meat has a high correlation with color parameters (L*—lightness, a*—redness, b*—yellowness and Chroma), as the pHfinal can affect the reaction of myoglobin to oxymyoglobin. The b* index in meat, on the other hand, may be related to the concentration of fat and/or the presence of carotenoids in the diet which can be affected by forage preservation processes, such as silage and hay, which significantly reduces by up to 80% carotenoids levels13. It is believed that the carotenoid concentrations in the diet of this study were similar between treatments and consequently in values of b* of meat. Values of a* and Chroma directly depend on the content and state of the heme pigments in the muscle, due to the chemical state of iron (Fe), playing an important role in meat color10. These parameters showed no significant difference between treatments, however, higher values of a* and Chroma in meat are desired, as a result of the increase in oxymyoglobin and decrease in metmyoglobin that provides the meat’s “bloom”. According to Dawson et al.17, the minimum critical value for meat luminosity (L*) is 34. Lower values of L are related to elevating pHfinal, which results in the high concentration of metmyoglobin, making the meat darker, which causes rejection by consumers for associating dark meat to as old meat.The meat’s presentation and more precisely its color is an important factor that can influence a consumer’s purchase decision, as it gives us the idea of freshness and meat’ quality. The L* and a* color parameters are the most representative for these characteristics18. Although in our research it did not have a significant effect on the color parameters, we can indicate that the meat obtained in this research would be well accepted by consumers, because Hopkins19 suggests that consumers will consider meat color acceptable when the L* value is equal to or exceeds 34, and a* value below 19 or equal to or exceeds 9.5 according to Khliji et al.18. In the present study, all values for L* remained above this aforementioned threshold and the values of a* remained within these values which suggests that meats from all diets and water supply levels had an acceptable color for consumers.When evaluating the chemical composition of meat, no significant differences were observed between treatments, except for the ash content, that remained above the average values found in the literature, which is 0.99–1.10%16. It is believed that because cactus pear is a rich source of Ca, Mg, K and with increasing level of cactus pear silage in the diet31, these minerals were consumed in larger amounts, which could have resulted in a higher proportion of minerals in the meat of animals that received 42% cactus pear silage.The lipid fatty acid profile in meat has a major impact on sensory properties and nutritional quality, influencing acceptance and health for consumers20,21. Intermittent water supply, cactus pear silage, and interaction between water supply and cactus pear silage did not influence most fatty acids present in the Longissimus lumborum muscle of the animals under study, except only a few saturated fatty acids e.g. docosanoic acid (C22:0), tricosanoic acid (C23:0), BCFA, anteiso-tridecanoic acid (C13:0 anteiso) and anteiso-pentadecanoic acid (C15:0 anteiso).Biohydrogenation of ruminal bacteria results in a circumstantial variety of fatty acids (FA), which will be absorbed in the intestine and later incorporated into the meat of goats. In addition to the diet and the biohydrogenation, the meat lipid profile can vary due to de novo synthesis, desaturation, duration of the feeding period and differences in pathways of various FA by the animal organism22.A high concentration of saturated fatty acids present in meat is not desirable, as there is evidence that saturated fatty acids, mainly C16:0, as well as myristic (C14:0) and lauric (C12:0) increase the blood cholesterol and low-density lipoproteins (LDL) concentration, due to interferences with hepatic LDL receptors23, however, in the studied treatments, there were no significant differences for these fatty acids. On the other hand, C18:0 has no impact on cholesterol levels, due to being poorly digested and easily desaturated to C18:1 by Δ9-desaturase24, present in the cell endoplasmic reticulum. This fatty acid is not harmful to health and is considered the only desirable SFA. As the levels of C18:0 in diets tend to be minimal, their main origin is the biohydrogenation of PUFA and de novo syntheses in diets with a high energy pattern25.In addition to carrying out the biohydrogenation process, ruminal bacteria synthesize a series of FA, mainly those of odd and branched chain, that comprise mainly the lipids of the bacterial membrane26,27, to maintain membrane fluidity. Linear odd-chains fatty acids are formed when propionyl-CoA, instead of acetyl-CoA, is used as a de novo synthesis initiator25. On the other hand, iso and anteiso FA are synthesized by the precursors branched-chain amino acids (valine, leucine, and isoleucine) and their corresponding branched- short-chain carboxylic acids (isobutyric, isovaleric and 2-methyl butyric acids)28.There is an increasing interest to study odd-and branched-chain fatty acids (OBCFAs) from animal products, mainly in milk due to its higher concentration compared to meat. Researchers reported that several OBCFAs have potential health benefits in humans29 as improved gut health30 and presenting anti-cancer activity31, as well as improve the sensory characteristics of the meat, providing a greater sensation of tenderness and juiciness, because BCFA content are associated with a less consistent fat in meat from lambs due to its lower melting point and its chain structure32.The FAs profile in the ruminal bacteria is largely composed by OBCFAs (C15:0; anteiso C15:0; iso C15:0; C17:0; iso C17:0; C17:1 and anteiso C17:0) in the bacteria membrane lipids24. Thus, the higher concentration of OBCFAs might be the result of the difference in the rumen bacterial populations induced by variation in the dietary carbohydrate, that is, a higher concentration of cellulolytic bacteria in relation to amylolytic bacteria, due to the high neutral detergent fiber (NDF) content in the diet with 0% cactus forage silage. It is also known that amylolytic bacteria produce more linear odd chain and anteiso FAs than iso FAs, whereas cellulolytic bacteria produce more iso FAs28,32. As the Tifton 85 grass hay-based diet had the highest neutral detergent fiber corrected for ash and protein (NDFap) and starch content (highest % of ground corn), the meat of those animals had higher concentrations of anteiso C15:0 and anteiso C13:0 compared to animals fed diets with the inclusion of cactus pear silage, also influencing the total sum of branched chain fatty acids.Although levels of intermittent water supply have generated punctual changes in tricosanoic acid (C23:0) SFA, the same was not observed for MUFA and PUFA, due to changes in the rumen environment, promoted by water restrictions, which were not sufficient to circumstantially modify biohydrogenation, resulting in similarities in concentrations of unsaturated fatty acids in goat meat.The animals subjected to 24 h of intermittent water supply (IWS) presented the highest concentration of C23:0 in relation to other treatments, which is interesting because it is involved in the synthesis of ceramide and reduces the risk of diabetes in humans33.The cactus pear has high non-fibrous carbohydrate (NFC) content (mainly pectin), having 59.5% high and medium rumen degradation carbohydrates which provide a higher production rate and removal of short-chain fatty acids and changes in rumen bacterial populations34. The inclusion of CPS resulted in a higher passage rate of digesta, affected biohydrogenation, and resulted in the escape of intermediate fatty acids isomers that are absorbed in the small intestine. Consequently, there was changing composition of fatty acids in the muscle of these animals, with a significant effect being observed only in the cis-13 C18:1. Furthermore, diets with high proportions of cactus pear silage (CPS), such as 42% CPS diet, can decrease ruminal pH and affect the final stages of biohydrogenation, resulting in the escape of intermediate fatty acids isomers, that are absorbed in the small intestine, which can explain the similarity of the C20:1 in 42% CPS diet from the Tifton hay-based diet, with differences between goat meat from 21% CPS diet and Tifton hay-based diet.Oleic acid (c9-C18:1) was the MUFA with the highest participation in the lipid profile of goat meat, which is interesting because it has a hypocholesterolemic effect, being a desirable fatty acid (DFA) for not reducing the serum high density lipoproteins (HDL) levels and thus prevent cardiovascular disease by reducing LDL levels35. The high concentrations of c9-C18:1 in ruminant meat come from the food intake, the effect of biohydrogenation, and mainly of the high activity of Δ9-desaturase, necessary for animal biosynthesis through desaturation of C18:0 to c9-C18:127. This fatty acid in the lipid profile of red meat varies between 30 and 43%36, confirming that the meat in the present study had a good concentration of this fatty acid.Much of unsaturated fatty acids, which have 18 carbons or 16 carbons, are largely converted to C18:0 and C16:0 through biohydrogenation, and when this process is not 100% completed, in addition to the PUFA that pass through this process intact, some product intermediates are formed, reaching the duodenum and are absorbed by the animal, in which significant amounts of cis and trans-monounsaturated, such as vaccenic fatty acid (t11-C18:1), reach the duodenum and are absorbed, later composing the muscle tissue22.The literature indicates that the precursor of conjugated linoleic acid (CLA) in the meat of animals is trans vaccenic acid (t11-C18:1), so the enzyme ∆9-desaturase, besides acting in the conversion of stearic into oleic fatty acid, also converts the trans-vaccenic acid to its corresponding CLA isomer, c9t11-C18:236. This pathway is more expressive in the mammary gland, and as the concentration of vaccenic acid (t11-C18:1) was not different, the concentration of CLA was not affected by the supply of silage and intermittent water supply, in the same way, that there are also no differences in the activity of ∆9-desaturase. Nevertheless, it is worth noting that in the human adipose tissue there is also the presence of ∆9-desaturase, and therefore, increased intake of vaccenic fatty acid could have the same beneficial effects associated with the intake of CLA, where the dietary vaccenic fatty acid shows 19–30% conversion rate37.Tifton hay is a natural source of n-3 fatty acids, mainly C18:3 n-3 with up to 20% participation in the lipid profile2, allowing a certain part of these PUFAs to be absorbed and increased in the tissue muscle, with 10 to 30% PUFAs in the diet generally escaping from biohydrogenation.Linoleic fatty acid (c9c12 C18:2) and α-linolenic acid (C18:3 n-3) are essential fatty acids for humans, that serve as precursors of the n-3 and n-6 pathways, distinct families, but synthesized by some of the same enzymes (∆4-desaturase, ∆5-desaturase, and ∆6-desaturase)25. Arachidonic fatty acid (C20:4 n-6) comes from elongation and desaturation of linoleic acid, where its concentrations, even close to that of its precursor, may indicate that there was a high activity of ∆6-desaturase (desaturation to γ-linolenic), elongase (elongation of γ-linolenic to dihomo-gamma-linolenic) and ∆5-desaturase. This fatty acid was influenced by the diets, presenting lower concentrations in the meat of animals fed the 42% cactus pear silage when compared to the Tifton hay diet (0% cactus pear silage).A higher concentration of long-chain PUFA n-3, docosahexaenoic (C22:6 n-3), was observed in the muscle of animals fed on Tifton hay. This was probably due to the high concentration of C18:3 n-3, precursor of C22:6 n-3, that the hay presents in relation to the cactus pear silage.The ratios and proportions of fatty acids are used to determine nutritional and nutraceutical values of the product or diet, and mainly, to indicate the cholesterolemic potential4. It is interesting that the n-6/n-3 ratio is low due to the pro-inflammatory properties of n-6; it is recommended to decrease its intake to assist in disease prevention38, while n-3 fatty acids are anti-inflammatory, antithrombotic, antiarrhythmic and reduce blood lipids, with vasodilating properties, being interesting that they present a higher proportion24. n-6 fatty acids tend to have a higher percentage in meat, and this directly influences the formation of n-3 isomers, since linoleic acid, when in excess, can reduce the synthesis of linolenic acid metabolites. The percentage of FA in one group can interfere with the metabolism of the other, reducing its incorporation into tissue lipids and altering its general biological effects38. Therefore, it is not recommended that the n-6/n-3 ratio be kept above 5 or 639, demonstrating that the averages of the current research remained acceptable.In relation to atherogenicity index (AI) and thrombogenicity index (TI), Ulbricht and Southgate39 proposed that sheep meat should have values of up to 1.0 and 1.58, respectively, and the lower the values for these indices in the lipid fraction, the greater the prevention of early stages of cardiovascular diseases. In the present study, the general averages observed were 0.29 for the AI, and 0.81 for the TI, although there were no significant differences, all treatments are within the recommended range, despite having been used as comparative standard to sheep, due to the absence of the proposed standard for goat meat.The h:H ratio did not differ for diets and water supply levels, but had an average of 1.90, below the reference value for meat products, which is 2.0. Values above 2.0 are recommended and favorable40, as it indicates a higher proportion of hypocholesterolemic fatty acids, that are beneficial to human health.The ∆9-desaturase enzyme that acts on both the mammary gland and adipose tissue, responsible for the transformation of SFA into unsaturated fatty acids (UFA), as well as in the endogenous conversion of CLA37 did not differ between treatments. On the other hand, the elongase showed less activity. Probably there was a greater “de novo” synthesis which resulted in a greater accumulation of palmitic fatty acid, and a reduction in the activity of the elongase enzyme.The crossbred goats demonstrated to present efficient mechanisms for adapting to water restrictions, especially when receiving feed with higher water content, such as cactus pear silage, being able to replace Tifton hay with 42% cactus pear silage in the diet for goats in confinement without negatively affecting the carcass traits and meat quality. Because, although these animals have shown some differences in the indices of tenderness and juiciness of their meats, however, all presented values of juiciness and tenderness compatible with meat extremely appreciated by the consumer market, and even goat meat showing some fatty acids with different concentrations induced by the supply of silage and water intermittence, the final lipid profile was appropriate to the health of consumers, observed by the absence of differences in the total concentrations of PUFA and in the main nutraceutical parameters (DFA, n-6/n-3; h:H; AI and TI).These results are relevant, indicating that goat feedlots in regions with low water availability may adopt strategies of lesser demand for drinking water and considerable concentrations of cactus pear silage in the diet, can reduce production costs without considerably affecting the product to be marketed, and therefore, provide higher profitability of the system. More