Ecology

Finkelman, S., Navarro, S., Rindner, M. & Dias, R. Effect of low pressure on the survival of Trogoderma granarium Everts, Lasioderma serricorne (F.) and Oryzaephilus surinamensis (L.) at 30°C. J. Stored. Prod. Res. 42, 23–30 (2006).Article 

Google Scholar 
Hosseininaveh, V., Bandani, A., Azmayeshfard, P., Hosseinkhani, S. & Kazzazi, M. Digestive proteolytic and amylolytic activities in Trogoderma granarium Everts (Dermestidae: Coleoptera). J. Stored. Prod. Res. 43, 515–522 (2007).Article 
CAS 

Google Scholar 
Burges, H. D. Development of the khapra beetle, Trogoderma granarium, in the lower part of its temperature range. J. Stored. Prod. Res. 44, 32–35 (2008).Article 

Google Scholar 
Hagstrum, D. W. & Subramanyam, B. Stored-Product Insect Resource 1–518 (AACC International Inc, 2009).Book 

Google Scholar 
Beal, R. S. Synopsis of the economic species of Trogoderma occurring in the United States with description of a new species (Coleoptera: Dermestidae). Ann. Entomol. Soc. Am. 49, 559–566 (1956).Article 

Google Scholar 
Day, C. & White, B. Khapra beetle, Trogoderma granarium interceptions and eradications in Australia and around the world. Crawley, School of Agricultural and Resource Economics, University of Western Australia, SARE Working paper 1609, (2016).Kerr, J. A. Khapra beetle returns. Pest Control 49, 24–25 (1981).
Google Scholar 
Stibick, J.N. New pest response guidelines: khapra beetle. US Department of Agriculture, Marketing and Regulatory Programs, Animal and Plant Health Inspection Service, Riverdale, pp. 114 (2009).Myers, S. W. & Hagstrum, D. W. Quarantine. In Stored Product Protection (eds Hagstrum, D. W. et al.) 297–304 (Kansas State University Agricultural Experiment Station and Cooperative Extension Service, 2012).
Google Scholar 
Athanassiou, C. G., Phillips, T. W. & Wakil, W. Biology and control of the khapra beetle, Trogoderma granarium, a major quarantine threat to global food security. Annu. Rev. Entomol. 64, 131–148 (2019).Article 
CAS 
PubMed 

Google Scholar 
Barak, A. V. Development of a new trap to detect and monitor khapra beetle (Coleoptera: Dermestidae). J. Econ. Entom. 82, 1470–1477 (1989).Article 

Google Scholar 
Gerken, A. R. & Campbell, J. F. Life history changes in Trogoderma variabile and T. inclusum due to mating delay with implications for mating disruption as a management tactic. Ecol. Evol. 8, 2428–2439 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Partida, G. J. & Strong, R. G. Comparative studies on the biologies of six species of Trogoderma: T variabile. Ann. Entomol. Soc. Am. 68, 115–125 (1975).Article 

Google Scholar 
Strong, R. G. Comparative studies on the biologies of six species of Trogoderma: T inclusum. Ann. Entomol. Soc. Am. 68, 91–104 (1975).Article 

Google Scholar 
Phillips, T. W., Pfannenstiel, L. & Hagstrum, D. Survey of Trogoderma species (Coleoptera: Dermestidae) associated with international trade of dried distiller’s grains and solubles in the USA. Julius Kühn Archiv. 463, 233–238 (2008).
Google Scholar 
Hadaway, A. The biology of the beetles, Trogoderma granarium Everts and Trogoderma versicolor (Creutz). Bull. Entomol. Res. 46, 781–796 (1956).Article 
CAS 

Google Scholar 
Phillips, T.W., Pfannenstiel, L. & Hagstrum, D. Survey of Trogoderma species (Coleoptera: Dermestidae) associated with international trade of dried distiller’s grains and solubles in the USA. In: Adler CS, Opit G, Fürstenau B, Müller-Blenkle C, Kern P, Arthur FH et al., editors. Proceedings of the 12th International Working Conference on Stored Product Protection; Vol. 1, Quedlinburg, Julius-Kühn-Archiv, pp. 233–238 (2018).Gorham, J.R. Insect and Mite Pests in Food: An Illustrated Key. Vol. 1 and 2. US Department of Agriculture, Agricultural Research Service (1991).Olson, R. L. O., Farris, R. E., Barr, N. B. & Cognato, A. I. Molecular identification of Trogoderma granarium (Coleoptera: Dermestidae) using the 16S gene. J. Pest Sci. 87, 701–710 (2014).Article 

Google Scholar 
Furui, S., Miyanoshita, A., Imamura, T., Minegishi, Y. & Kokutani, R. Qualitative real-time PCR identification of the khapra beetle, Trogoderma granarium (Coleoptera: Dermestidae). Appl. Entomol. Zool. 54, 101–107 (2019).Article 
CAS 

Google Scholar 
Rako, L. et al. A LAMP (loop-mediated isothermal amplification) test for rapid identification of Khapra beetle (Trogoderma granarium). Pest Manag. Sci. 77, 5509–5521 (2021).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Castañé, C., Agustí, N., del Estal, P. & Riudavets, J. Survey of Trogoderma spp. in Spanish mills and warehouses. J. Stored. Prod. Res. 88, 101661 (2020).Article 

Google Scholar 
Trujillo-González, et al. Detection of khapra beetle environmental DNA using portable technologies in Australian biosecurity. Front. Insect Sci. 2, e795379 (2022).Article 

Google Scholar 
Svec, D., Tichopad, A., Novosadova, V., Pfaffl, M. W. & Kubista, M. How good is a PCR efficiency estimate: Recommendations for precise and robust qPCR efficiency assessments. Biomol. Detect. Quantif. 3, 9–16 (2015).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Taylor, S. C. et al. The Ultimate qPCR experiment: Producing publication quality, reproducible data the first time. Trends Biotechnol. 37, 761–774 (2019).Article 
CAS 
PubMed 

Google Scholar 
Van Holm, W. et al. A viability quantitative PCR dilemma: Are longer amplicons better?. Appl. Environ. Microbiol. 87, e0265320 (2021).Article 
PubMed 

Google Scholar 
Ratnasingham, S. & Hebert, P. D. N. BOLD: The barcode of life data system (wwwbarcodinglifeorg). Mol. Ecol. Notes 7, 355–364 (2007).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Wittwer, C. T. & Kusakawa, N. Real-time PCR. In Molecular microbiology: Diagnostic principles and practice (eds Persing, D. H. et al.) 71–84 (ASM Press, 2004).
Google Scholar 
Stewart, D. et al. A needle in a haystack: A multigene TaqMan assay for the detection of Asian gypsy moths in bulk pheromone trap samples. Biol. Invasions 21, 1843–1856 (2019).Article 

Google Scholar 
Butterwort, V. et al. A DNA extraction method for insects from sticky traps: Targeting a low abundance pest, Phthorimaea absoluta (Lepidoptera: Gelechiidae), in mixed species communities. J. Econ. Entom. 115, 844–851 (2022).Article 
CAS 
PubMed 

Google Scholar 
Carew, M. E., Coleman, R. A. & Hoffmann, A. A. Can non-destructive DNA extraction of bulk invertebrate samples be used for metabarcoding?. PeerJ 6, e4980 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Domingue, M.J. et al. Outcome of direct competition between Trogoderma granarium and Trogoderma inclusum over varying commodities, temperatures, and experimental duration. In Submission to Scientific Reports.Zieritz, A. et al. Development and evaluation of hotshot protocols for cost- and time-effective extraction of PCR-ready DNA from single freshwater mussel larvae (Bivalvia: Unionida). J. Molluscan Stud. 84, 198–201 (2018).Article 

Google Scholar 
Djoumad, A. et al. Development of a qPCR-based method for counting overwintering spruce budworm (Choristoneura fumiferana) larvae collected during fall surveys and for assessing their natural enemy load: A proof-of-concept study. Pest Manag. Sci. 78, 336–343 (2022).Article 
CAS 
PubMed 

Google Scholar 
Chen, H., Rangasamy, M., Tan, S. Y., Wang, H. & Siegfried, B. D. Evaluation of five methods for total DNA extraction from western corn rootworm beetles. PLoS ONE 5, e11963 (2010).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Beckmann, J. S. & Soller, M. Restriction fragment length polymorphisms in genetic improvement: Methodologies, mapping and costs. Theor. Appl. Genet. 67, 35–43 (1983).Article 
CAS 
PubMed 

Google Scholar 
Arimoto, M., Satoh, M., Uesugi, R. & Osakabe, M. PCR-RFLP analysis for identification of Tetranychus spider mite species (Acari: Tetranychidae). J. Econ. Entom. 106, 661–668 (2013).Article 
CAS 
PubMed 

Google Scholar 
Vezenegho, S. B. et al. Discrimination of 15 Amazonian anopheline mosquito species by polymerase chain reaction—Restriction fragment length polymorphism. J. Med. Entomol. 59, 1060–1064 (2022).Article 
CAS 
PubMed 

Google Scholar 
Beal, R. S. Annotated checklist of Nearctic Dermestidae with revised key to the genera. Coleopt. Bull. 57, 391–404 (2003).Article 

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

Google Scholar 
Liu, H. & Mottern, J. An old remedy for a new problem? Identification of Ooencyrtus kuvanae (Hymenoptera: Encyrtidae), an egg parasitoid of Lycorma delicatula (Hemiptera: Fulgoridae) in North America. J. Insect Sci. 17, 18 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Simon, C. et al. Evolution, weighing, and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Ann. Entomol. Soc. Am. 87, 651–701 (1994).Article 
CAS 

Google Scholar 
Dowton, M. & Austin, A. D. Evidence for AT-transversion bias in wasp (Hymenoptera: Symphyta) mitochondrial genes and its implications for the origin of parasitism. J. Mol. Evol. 44, 398–405 (1997).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Untergasser, A. et al. Primer3—New capabilities and interfaces. Nucleic Acids Res. 40, e115 (2012).Article 
CAS 
PubMed 
PubMed Central 

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

Google Scholar 
Süss, B., Flekna, G., Wagner, M. & Hein, I. Studying the effect of single mismatches in primer and probe binding regions on amplification curves and quantification in real-time PCR. J. Microbiol. Methods 76, 316–319 (2009).Article 
PubMed 

Google Scholar 
Stadhouders, R. et al. The effect of primer-template mismatches on the detection and quantification of nucleic acids using the 5’ nuclease assay. J. Mol. Diagn. 12, 109–117 (2010).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Stewart, D. et al. A multi-species TaqMan PCR assay for the identification of Asian gypsy moths (Lymantria spp.) and other invasive Lymantriines of biosecurity concern to North America. PLoS ONE 11, e0160878 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Bustin, S. A. et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 55, 1–12 (2009).Article 

Google Scholar  More

Breed differences in animal feed conversion and economic trait performanceThroughout the feed intake measurement period, summary statistics shows animals on test had an average DMI of 1.11 kg/d (SD = 0.18), ADG of 0.27 kg/d (SD = 0.1), FCR of 4.04 kg of DMI/ Kg of ADG (SD = 0.1), start weight of 29.60 kg (SD = 3.7), final live weight of 46.00 kg (SD = 2.9), carcass weight of 20.20 kg (SD = 1.6), and a KO% of 44.1% (SD = 2.3). Average daily gain (P = 0.005), FCR (P = 0.035), CW (P  More

Semin, G. R., Scandurra, A., Baragli, P., Lanatà, A. & D’Aniello, B. Inter- and intra-species communication of emotion: Chemosignals as the neglected medium. Animals 9, 887 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Désiré, L., Boissy, A. & Veissier, I. Emotions in farm animals: A new approach to animal welfare in applied ethology. Behav. Processes 60, 165–180 (2002).Article 
PubMed 

Google Scholar 
Briefer, E. F. & Le Comber, S. Vocal expression of emotions in mammals: mechanisms of production and evidence. J. Zool. 288, 1–20 (2012).Article 

Google Scholar 
Jardat, P. & Lansade, L. Cognition and the human–animal relationship: a review of the sociocognitive skills of domestic mammals toward humans. Anim. Cogn. 25, 369–384 (2022).Article 
PubMed 

Google Scholar 
Galvan, M. & Vonk, J. Man’s other best friend: domestic cats (F. silvestris catus) and their discrimination of human emotion cues. Anim. Cogn. 19, 193–205 (2016).Nawroth, C., Albuquerque, N., Savalli, C., Single, M. S. & McElligott, A. G. Goats prefer positive human emotional facial expressions. R. Soc. Open Sci. 5, 180491 (2018).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Proops, L., Grounds, K., Smith, A. V. & McComb, K. Animals remember previous facial expressions that specific humans have exhibited. Curr. Biol. 28, 1428-1432.e4 (2018).Article 
CAS 
PubMed 

Google Scholar 
Smith, A. V., Proops, L., Grounds, K., Wathan, J. & McComb, K. Functionally relevant responses to human facial expressions of emotion in the domestic horse (Equus caballus). Biol. Lett. 12, 20150907 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Siniscalchi, M., D’Ingeo, S., Fornelli, S. & Quaranta, A. Lateralized behavior and cardiac activity of dogs in response to human emotional vocalizations. Sci. Rep. 8, 77 (2018).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Siniscalchi, M., D’Ingeo, S. & Quaranta, A. Orienting asymmetries and physiological reactivity in dogs’ response to human emotional faces. Learn. Behav. 46, 574–585 (2018).Article 
PubMed 

Google Scholar 
Smith, A. V. et al. Domestic horses (Equus caballus) discriminate between negative and positive human nonverbal vocalisations. Sci. Rep. 8, 13052 (2018).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Trösch, M. et al. Horses categorize human emotions cross-modally based on facial expression and non-verbal vocalizations. Animals 9, 862 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Quaranta, A., D’ingeo, S., Amoruso, R. & Siniscalchi, M. Emotion recognition in cats. Animals 10, 1107 (2020).Nakamura, K., Takimoto-Inose, A. & Hasegawa, T. Cross-modal perception of human emotion in domestic horses (Equus caballus). Sci. Rep. 8, 8660 (2018).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Albuquerque, N. et al. Dogs recognize dog and human emotions. Biol. Lett. 12, 20150883 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Briefer, E. F. Vocal contagion of emotions in non-human animals. Proc. R. Soc. B Biol. Sci. 285 (2018).Sabiniewicz, A., Tarnowska, K., Świątek, R., Sorokowski, P. & Laska, M. Olfactory-based interspecific recognition of human emotions: Horses (Equus ferus caballus) can recognize fear and happiness body odour from humans (Homo sapiens). Appl. Anim. Behav. Sci. 230, 105072 (2020).Article 

Google Scholar 
Baba, C., Kawai, M. & Takimoto-Inose, A. Are horses (Equus caballus) sensitive to human emotional cues?. Animals 9, 630 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
Brennan, P. A. & Kendrick, K. M. Mammalian social odours: attraction and individual recognition. Philos. Trans. R. Soc. B Biol. Sci. 361, 2061–2078 (2006).Saslow, C. A. Understanding the perceptual world of horses. Appl. Anim. Behav. Sci. 78, 209–224 (2002).Article 

Google Scholar 
Péron, F., Ward, R. & Burman, O. Horses (Equus caballus) discriminate body odour cues from conspecifics. Anim. Cogn. 17, 1007–1011 (2014).Article 
PubMed 

Google Scholar 
Krueger, K. & Flauger, B. Olfactory recognition of individual competitors by means of faeces in horse (Equus caballus). Anim. Cogn. 14, 245–257 (2011).Article 
PubMed 

Google Scholar 
Boissy, A., Terlouw, C. & Le Neindre, P. Presence of cues from stressed conspecifics increases reactivity to aversive events in cattle: evidence for the existence of alarm substances in urine. Physiol. Behav. 63, 489–495 (1998).Article 
CAS 
PubMed 

Google Scholar 
Vieuille-Thomas, C. & Signoret, J. P. Pheromonal transmission of an aversive experience in domestic pig. J. Chem. Ecol. 18, 1551–1557 (1992).Article 
CAS 
PubMed 

Google Scholar 
Siniscalchi, M., D’Ingeo, S. & Quaranta, A. The dog nose ‘KNOWS’ fear: Asymmetric nostril use during sniffing at canine and human emotional stimuli. Behav. Brain Res. 304, 34–41 (2016).Article 
PubMed 

Google Scholar 
Calvi, E. et al. The scent of emotions: A systematic review of human intra- and interspecific chemical communication of emotions. Brain Behav. 10 (2020).de Groot, J. H. B., Semin, G. R. & Smeets, M. A. M. On the communicative function of body odors: A theoretical integration and review. Perspect. Psychol. Sci. 12, 306–324 (2017).Article 
PubMed 

Google Scholar 
de Groot, J. H. B. et al. A sniff of happiness. Psychol. Sci. 26, 684–700 (2015).Article 
PubMed 

Google Scholar 
Destrez, A. et al. Male mice and cows perceive human emotional chemosignals: A preliminary study. Anim. Cogn. 24, 1205–1214 (2021).Article 
PubMed 

Google Scholar 
Wilson, Id. Dogs can discriminate between human baseline and psychological stress condition odours. PLoS ONE 17, e0274143 (2022).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
D’Aniello, B., Semin, G. R., Alterisio, A., Aria, M. & Scandurra, A. Interspecies transmission of emotional information via chemosignals: from humans to dogs (Canis lupus familiaris). Anim. Cogn. 21, 67–78 (2018).Article 
PubMed 

Google Scholar 
D’Aniello, B. et al. Sex differences in the behavioral responses of dogs exposed to human chemosignals of fear and happiness. Anim. Cogn. 24, 299–309 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Siniscalchi, M., d’Ingeo, S., Quaranta, A., D’Ingeo, S. & Quaranta, A. Lateralized emotional functioning in domestic animals. Appl. Anim. Behav. Sci. 237, 105282 (2021).Article 

Google Scholar 
Rogers, L. & Vallortigara, G. Lateralized Brain Functions: Methods in Human and Non-Human Species. vol. 122 (2017).D’Ingeo, S. et al. Horses associate individual human voices with the valence of past interactions: A behavioural and electrophysiological study. Sci. Rep. 9, 11568 (2019).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Siniscalchi, M., Padalino, B., Aubé, L. & Quaranta, A. Right-nostril use during sniffing at arousing stimuli produces higher cardiac activity in jumper horses. Laterality 20, (2015).De Boyer Des Roches, A., Richard-Yris, M.-A. A., Henry, S., Ezzaouïa, M. & Hausberger, M. Laterality and emotions: Visual laterality in the domestic horse (Equus caballus) differs with objects’ emotional value. Physiol. Behav. 94, 487–490 (2008).Albrecht, J. et al. Smelling chemosensory signals of males in anxious versus nonanxious condition increases state anxiety of female subjects. Chem. Senses 36, 19–27 (2011).Article 
CAS 
PubMed 

Google Scholar 
Derrickson, S. Sinister (VF)—YouTube. (Wild Bunch SA, 2001).Friard, O. & Gamba, M. BORIS: a free, versatile open-source event-logging software for video/audio coding and live observations. Methods Ecol. Evol. 7, 1325–1330 (2016).Article 

Google Scholar 
R Core Team. R: A language and environment for statistical computing (2021).Wickham, H. Ggplot2: Elegant graphics for data analysis. (2016).Hothorn, T., Winell, H., Hornik, K., van de Wiel, M. A. & Zeileis, A. coin: Conditional inference procedures in a permutation test framework (2021).Brooks, M. E. et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9, 378–400 (2017).Article 

Google Scholar 
Hartig, F. DHARMa: Residual diagnostics for hierarchical (multi-level/mixed) regression models (2021).Lenth, R. V. emmeans: Estimated Marginal Means, aka Least-Squares Means. (2022).Hothersall, B., Harris, P., Sörtoft, L. & Nicol, C. J. Discrimination between conspecific odour samples in the horse (Equus caballus). Appl. Anim. Behav. Sci. 126, 37–44 (2010).Article 

Google Scholar 
Smeets, M. A. M. et al. Chemical fingerprints of emotional body odor. Metabolites 10, (2020).Sabiniewicz, A. et al. A preliminary investigation of interspecific chemosensory communication of emotions: Can Humans (Homo sapiens) recognise fear- and non-fear body odour from horses (Equus ferus caballus). Animal 11, 3499 (2021).Article 

Google Scholar 
Zhou, W. & Chen, D. Entangled chemosensory emotion and identity: Familiarity enhances detection of chemosensorily encoded emotion. Soc. Neurosci. 6, 270–276 (2011).Article 
PubMed 

Google Scholar 
Starling, M., McLean, A. & McGreevy, P. The contribution of equitation science to minimising horse-related risks to humans. Animal 6, 15 (2016).Article 

Google Scholar 
Basile, M. et al. Socially dependent auditory laterality in domestic horses (Equus caballus). Anim. Cogn. 12, 611–619 (2009).Article 
PubMed 

Google Scholar 
Hatfield, E., Cacioppo, J. T. & Rapson, R. L. Emotional contagion. Curr. Dir. Psychol. Sci. 2, 96–100 (1993).Article 

Google Scholar 
Austin, N. P. & Rogers, L. J. Limb preferences and lateralization of aggression, reactivity and vigilance in feral horses Equus caballus. Anim. Behav. 83, 239–247 (2012).Article 

Google Scholar 
Larose, C., Richard-Yris, M.-A., Hausberger, M. & Rogers, L. J. Laterality of horses associated with emotionality in novel situations Laterality Asymmetries Body. Brain Cogn. 11, 355–367 (2006).
Google Scholar 
Farmer, K., Krüger, K., Byrne, R. W. & Marr, I. Sensory laterality in affiliative interactions in domestic horses and ponies (Equus caballus). Anim. Cogn. 21, 631–637 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Chen, D. & Haviland-Jones, J. Human olfactory communication of emotion. Percept. Mot. Skills 91, 771–781 (2000).Article 
CAS 
PubMed 

Google Scholar 
de Groot, J. H. B., Semin, G. R. & Smeets, M. A. M. Chemical communication of fear: A case of male-female asymmetry. J. Exp. Psychol. Gen. 143, 1515–1525 (2014).Article 
PubMed 

Google Scholar 
Marinier, S. L., Alexander, A. J. & Waring, G. H. Flehmen behaviour in the domestic horse: Discrimination of conspecific odours. Appl. Anim. Behav. Sci. 19, 227–237 (1988).Article 

Google Scholar 
Lansade, L., Pichard, G. & Leconte, M. Sensory sensitivities: Components of a horse’s temperament dimension. Appl. Anim. Behav. Sci. 114, 534–553 (2008).Article 

Google Scholar 
Lansade, L., Bouissou, M. F. & Erhard, H. W. Fearfulness in horses: A temperament trait stable across time and situations. Appl. Anim. Behav. Sci. 115, 182–200 (2008).Article 

Google Scholar 
Hoenen, M., Wolf, O. T. & Pause, B. M. The impact of stress on odor perception. https://doi.org/10.1177/030100661668870746,366-376 (2017).Article 
PubMed 

Google Scholar 
Rørvang, M. V., Nicova, K. & Yngvesson, J. Horse odor exploration behavior is influenced by pregnancy and age. Front. Behav. Neurosci. 16, 295 (2022).Article 

Google Scholar 
Doty, R. L. & Cameron, E. L. Sex differences and reproductive hormone influences on human odor perception. Physiol. Behav. 97, 213–228 (2009).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar  More

Study populations and sequencing strategyDNA libraries were prepared for 1261 D. sylvestris individuals from 115 populations (5–20 individuals per population) under a modified protocol49 of the Illumina Nextera DNA library preparation kit (Supplementary Methods S1.1, Supplementary Data 1). Individuals were indexed with unique dual-indexes (IDT Illumina Nextera 10nt UDI – 384 set) from Integrated DNA Technologies Co, to avoid index-hopping50. Libraries were sequenced (150 bp paired-end sequencing) in four lanes of an Illumina NovaSeq 6000 machine at Novogene Co. This resulted in an average coverage of ca. 2x per individual. Sequenced individuals were trimmed for adapter sequences (Trimmomatic version 0.3551), mapped (BWA-MEM version 0.7.1752,53) against a reference assembly54 (ca. 440 Mb), had duplicates marked and removed (Picard Toolkit version 2.0.1; http://broadinstitute.github.io/picard), locally realigned around indels (GATK version 3.555), recalibrated for base quality scores (ATLAS version 0.956) and had overlapping read pairs clipped (bamUtil version 1.0.1457) (Supplementary Methods S1.1). Population genetic analyses were performed on the resultant BAM files via genotype likelihoods (ANGSD version 0.93358 and ATLAS versions 0.9–1.056), to accommodate the propagation of uncertainty from the raw sequence data to population genetic inference.Population genetic structure and biogeographic barriersTo investigate the genetic structure of our samples (Fig. 2A, Supplementary Fig. S2), we performed principal component analyses (PCA) on all 1261 samples (“full” dataset) via PCAngsd version 0.9859, following conversion of the mapped sequence data to ANGSD genotype likelihoods in Beagle format (Supplementary Methods S1.2). To visualise PCA results in space (Supplementary Fig. S4), individuals’ principal components were projected on a map, spatially interpolated (linear interpolation, akima R package version 0.6.260) and had the first two principal components represented as green and blue colour channels. Given that uneven sampling can bias the inference of structure in PCA, PCA was also performed on a balanced dataset comprising a common, down-sampled size of 125 individuals per geographic region (“balanced” dataset; Fig. 2B, Supplementary Fig. S3; Supplementary Methods S1.2; Supplementary Data 1). Individual admixture proportions and ancestral allele frequencies were estimated using PCAngsd (-admix model) for K = 2–6, using the balanced dataset to avoid potential biases related to imbalanced sampling22,23 and an automatic search for the optimal sparseness regularisation parameter (alpha) soft-capped to 10,000 (Supplementary Methods S1.2). To visualise ancestry proportions in space, population ancestry proportions were spatially interpolated (kriging) via code modified from Ref. 61 (Supplementary Fig. S5).To test if between-lineage admixture underlies admixture patterns inferred by PCAngsd or if the data is better explained by alternative scenarios such as recent bottlenecks, we used chromosome painting and patterns of allele sharing to construct painting palettes via the programmes MixPainter and badMIXTURE (unlinked model)28 and compared this to the PCAngsd-inferred palettes (Fig. 2B, C; Supplementary Methods S1.2). We referred to patterns of residuals between these palettes to inform of the most likely underlying demographic scenario. For assessing Alpine–Balkan palette residuals (and hence admixture), 65 individuals each from the French Alps (inferred as pure Alpine ancestry in PCAngsd), Monte Baldo (inferred with both Alpine and Balkan ancestries in PCAngsd) and Julian Alps (inferred as pure Balkan ancestry in PCAngsd) were analysed under K = 2 in PCAngsd and badMIXTURE (Fig. 2C). For assessing Apennine–Balkan admixture, 22 individuals each from the French pre-Alps (inferred as pure Apennine ancestry in PCAngsd), Tuscany (inferred with both Apennine and Balkan ancestries in PCAngsd) and Julian Alps (inferred as pure Balkan ancestry in PCAngsd) were analysed under K = 2 in PCAngsd and badMIXTURE.To construct a genetic distance tree (Supplementary Fig. S1), we first calculated pairwise genetic distances between 549 individuals (5 individuals per population for all populations) using ATLAS, employing a distance measure (weight) reflective of the number of alleles differing between the genotypes (Supplementary Methods S1.2; Supplementary Data 1). A tree was constructed from the resultant distance matrix via an initial topology defined by the BioNJ algorithm with subsequent topological moves performed via Subtree Pruning and Regrafting (SPR) in FastME version 2.1.6.162. This matrix of pairwise genetic distances was also used as input for analyses of effective migration and effective diversity surfaces in EEMS25. EEMS was run setting the number of modelled demes to 1000 (Fig. 2A, Supplementary Fig. S8). For each case, ten independent Markov chain Monte Carlo (MCMC) chains comprising 5 million iterations each were run, with a 1 million iteration burn-in, retaining every 10,000th iteration. Biogeographic barriers (Fig. 2A, Supplementary Fig. S7) were further identified via applying Monmonier’s algorithm24 on a valuated graph constructed via Delauney triangulation of population geographic coordinates, with edge values reflecting population pairwise FST; via the adegenet R package version 2.1.163. FST between all population pairs were calculated via ANGSD, employing a common sample size of 5 individuals per population (Supplementary Fig. S6; Supplementary Methods S1.2; Supplementary Data 1). 100 bootstrap runs were performed to generate a heatmap of genetic boundaries in space, from which a weighted mean line was drawn (Supplementary Fig. S7). All analyses in ANGSD were performed with the GATK (-GL 2) model, as we noticed irregularities in the site frequency spectra (SFS) with the SAMtools (-GL 1) model similar to that reported in Ref. 58 with particular BAM files. All analyses described above were performed on the full genome.Ancestral sequence reconstructionTo acquire ancestral states and polarise site-frequency spectra for use in the directionality index ψ and demographic inference, we reconstructed ancestral genome sequences at each node of the phylogenetic tree of 9 Dianthus species: D. carthusianorum, D. deltoides, D. glacialis, D. sylvestris (Apennine lineage), D. lusitanus, D. pungens, D. superbus alpestris, D. superbus superbus, and D. sylvestris (Alpine lineage). This tree topology was extracted from a detailed reconstruction of Dianthus phylogeny based on 30 taxa by Fior et al. (Fior, Luqman, Scharmann, Zemp, Zoller, Pålsson, Gargano, Wegmann & Widmer; paper in preparation) (Supplementary Methods S1.3). For ancestral sequence reconstruction, one individual per species was sequenced at medium coverage (ca. 10x), trimmed (Trimmomatic), mapped against the D. sylvestris reference assembly (BWA-MEM) and had overlapping read pairs clipped (bamUtil) (Supplementary Methods S1.3). For each species, we then generated a species-specific FASTA using GATK FastaAlternateReferenceMaker. This was achieved by replacing the reference bases at polymorphic sites with species-specific variants as identified by freebayes64 (version 1.3.1; default parameters), while masking (i.e., setting as “N”) sites (i) with zero depth and (ii) that didn’t pass the applied variant filtering criteria (i.e., that are not confidently called as polymorphic; Supplementary Methods S1.3). Species FASTA files were then combined into a multi-sample FASTA. Using this, we probabilistically reconstructed ancestral sequences at each node of the tree via PHAST (version 1.4) prequel65, using a tree model produced by PHAST phylofit under a REV substitution model and the specified tree topology (Supplementary Methods S1.3). Ancestral sequence FASTA files were then generated from the prequel results using a custom script.Expansion signalTo calculate the population pairwise directionality index ψ for the Alpine lineage, we utilised equation 1b from Peter and Slatkin (2013)31, which defines ψ in terms of the two-population site frequency spectrum (2D-SFS) (Supplementary Methods S1.4). 2D-SFS between all population pairs (10 individuals per population; Supplementary Data 1) were estimated via ANGSD and realSFS66 (Supplementary Methods S1.4), for unfolded spectra. Unfolding of spectra was achieved via polarisation with respect to the ancestral state of sites defined at the D. sylvestris (Apennine lineage) – D. sylvestris (Alpine lineage) ancestral node. Correlation of pairwise ψ and (great-circle) distance matrices was tested via a Mantel test (10,000 permutations). To infer the geographic origin of the expansion (Fig. 3), we employed a time difference of arrival (TDOA) algorithm following Peter and Slatkin (2013);31 performed via the rangeExpansion R package version 0.0.0.900031,67. We further estimated the strength of the founder of this expansion using the same package.Demographic inferenceTo evaluate the demographic history of D. sylvestris, a set of candidate demographic models was formulated. To constrain the topology of tested models, we first inferred the phylogenetic tree of the three identified evolutionary lineages of D. sylvestris (Alpine, Apennine and Balkan) as embedded within the larger phylogeny of the Eurasian Dianthus clade (note that the phylogeny from Fior et al. (Fior, Luqman, Scharmann, Zemp, Zoller, Pålsson, Gargano, Wegmann & Widmer; paper in preparation) excludes Balkan representatives of D. sylvestris). Trees were inferred based on low-coverage whole-genome sequence data of 1–2 representatives from each D. sylvestris lineage, together with whole-genome sequence data of 7 other Dianthus species, namely D. carthusianorum, D. deltoides, D. glacialis, D. lusitanus, D. pungens, D. superbus alpestris and D. superbus superbus, that were used to root the D. sylvestris clade (Supplementary Methods S1.5). We estimated distance-based phylogenies using ngsDist68 that accommodates genotype likelihoods in the estimation of genetic distances (Supplementary Methods S1.5). Genetic distances were calculated via two approaches: (i) genome-wide and (ii) along 10 kb windows. For the former, 110 bootstrap replicates were calculated by re-sampling over similar-sized genomic blocks. For the alternative strategy based on 10 kb windows, window trees were combined using ASTRAL-III version 5.6.369 to generate a genome-wide consensus tree accounting for potential gene tree discordance (Supplementary Methods S1.5). Trees were constructed from matrices of genetic distances from initial topologies defined by the BioNJ algorithm with subsequent topological moves performed via Subtree Pruning and Regrafting (SPR) in FastME version 2.1.6.162. We rooted all resultant phylogenetic trees with D. deltoides as the outgroup70. Both approaches recovered a topology with the Balkan lineage diverging prior to the Apennine and Alpine lineages (Supplementary Fig. S9). This taxon topology for D. sylvestris was supported by high ASTRAL-III posterior probabilities ( >99%), ASTRAL-III quartet scores ( >0.5) and bootstrap values ( >99%). Topologies deeper in the tree were less well-resolved (with quartet scores More

Modular components of the DFTe energy functionalThe central ingredient of DFTe is an energy functional E, assembled according to Eq. (1). The methodology of DFTe can be understood by inspecting the dispersal and environmental energies in Eqs. (2) and (3) without interactions. In our first case study, illustrated in Fig. 2 and Supplementary Fig. 2, we demonstrate that equation (3), in conjunction with Eq. (2), can realistically describe the influence of the environment on species’ distributions. Mechanisms that alter the trade-off between dispersal and environment can be introduced as part of Eint. For instance, back reactions on the environment could be modelled with a bifunctional Ebr[Venv, n] that yields the equilibrated modified environment ({V}_{s}^{{{{{{{{rm{env}}}}}}}}}+delta {E}_{{{{{{{{rm{br}}}}}}}}}[{{{{{{{{bf{V}}}}}}}}}^{{{{{{{{rm{env}}}}}}}}},{{{{{{{bf{n}}}}}}}}]/delta {n}_{s}({{{{{{{bf{r}}}}}}}})), cf. Eq. (5).In the following we make explicit the interaction and resource energies that enter Eq. (1) and are used in our case studies of Figs. 2–7. We let Eint[n] include all possible bipartite interactions$${E}_{gamma }[{{{{{{{bf{n}}}}}}}}]=mathop{sum }limits_{{s,{s}^{{prime} }!=!1}atop {{s}^{{prime} }ne s}}^{S}{int}_{A}({{{{{{{rm{d}}}}}}}}{{{{{{{bf{r}}}}}}}})({{{{{{{rm{d}}}}}}}}{{{{{{{{bf{r}}}}}}}}}^{{prime} }){n}_{s}{({{{{{{{bf{r}}}}}}}})}^{{alpha }_{s}},{gamma }_{s{s}^{{prime} }}({{{{{{{bf{r}}}}}}}},, {{{{{{{{bf{r}}}}}}}}}^{{prime} }){n}_{{s}^{{prime} }}{({{{{{{{{bf{r}}}}}}}}}^{{prime} })}^{{beta }_{{s}^{{prime} }}},$$
(6)
which include amensalism, commensalism, mutualism, and so forth. Here, ({alpha }_{s},, {beta }_{{s}^{{prime} }}ge 0), and the interaction kernels ({gamma }_{s{s}^{{prime} }}) are assembled from fitness proxies of species s and ({s}^{{prime} }) (Supplementary Table 1). Higher-order interactions can be introduced, for example, through (i) terms like ({n}_{s},{gamma }_{s{s}^{{prime} }},{n}_{{s}^{{prime} }},{gamma }_{{s}^{{prime} }{s}^{{primeprime} }}^{{prime} },{n}_{{s}^{{primeprime} }}) that build on pairwise interactions or (ii) genuinely multipartite expressions like ({gamma }_{s{s}^{{prime} }{s}^{{primeprime} }}{n}_{s},,{n}_{{s}^{{prime} }},{n}_{{s}^{{primeprime} }}). Multi-partite interactions based on bipartite interactions do not seem to be an uncommon scenario48. However, there may be systems where nonzero coefficients ({gamma }_{s{s}^{{prime} }{s}^{{primeprime} }}) couple all species. This poses a challenge for mechanistic theories in general. Then, ‘simpler subsystems’ that have to be included in the DFTe workflow of Fig. 1a can only refer to situations where other energy components are absent, such as resource terms or complex environments. For example, the coefficients ({gamma }_{s{s}^{{prime} }{s}^{{primeprime} }}) could be extracted in an experiment with a controlled simple environment and then used to model the interacting species in a real-world setting. For (({alpha }_{s},, {beta }_{{s}^{{prime} }})=(1,1)) we identify the contact interaction in physics as ({gamma }_{s{s}^{{prime} }}propto delta ({{{{{{{bf{r}}}}}}}}-{{{{{{{{bf{r}}}}}}}}}^{{prime} })) with the two-dimensional delta function δ( ), while the Coulomb interaction amounts to setting ({gamma }_{s{s}^{{prime} }}propto 1/|{{{{{{{bf{r}}}}}}}}-{{{{{{{{bf{r}}}}}}}}}^{{prime} }|). The mechanistic effect of these interaction kernels on the density distributions is the same in ecology as it is in physics—a mathematical insight that inspired us to build ecological analogues to the phenomenology of quantum gases, which feature functionals of the kind in Eq. (6). Note that we do not introduce any quantum effects into ecology despite the fact that the mathematical structure of DFTe is borrowed in part from quantum physics. While the contact interaction is a suitable candidate for plants and especially microbes52, we expect long-range interactions (for example, repulsion of Coulomb type) to be more appropriate for species with long-range sensors, such as eyes. Both types of interactions feature in describing the ecosystems addressed in this work.In a natural setting the equilibrium abundances are ultimately constrained by the accessible resources. It is within these limits of resource availability that environment as well as intra- and inter-specific interactions can shape the density distributions. An energy term for penalising over- and underconsumption of resources is thus of central importance. Each species consumes resources from some of the K provided resources, indexed by k. A subset of species consumes the locally available resource density ρk(r) according to the resource requirements νks, which represent the absolute amount of resource k consumed by one individual (or aggregated constituent) of species s. The simple quadratic functional$${E}_{{{{{{{{rm{Res}}}}}}}}}[{{{{{{{bf{n}}}}}}}}]={int}_{A}({{{{{{{rm{d}}}}}}}}{{{{{{{bf{r}}}}}}}})mathop{sum }limits_{k=1}^{K}{{{{{{{{mathcal{L}}}}}}}}}_{k}left({{{{{{{bf{n}}}}}}}},, {rho }_{k}right)equiv zeta {int}_{A}({{{{{{{rm{d}}}}}}}}{{{{{{{bf{r}}}}}}}})mathop{sum }limits_{k=1}^{K}{w}_{k}({{{{{{{bf{r}}}}}}}}){left[mathop{sum }limits_{s=1}^{S}{nu }_{ks}{n}_{s}({{{{{{{bf{r}}}}}}}})-{rho }_{k}({{{{{{{bf{r}}}}}}}})right]}^{2}$$
(7)
proves appropriate. Here, νksns is the portion of resource density ρk that is consumed by species s. That is, νks  > 0 indicates that s requires resource k. If Eq. (7) is the total energy functional, then a single-species system with a single resource equilibrates with density n1(r) = ρ1(r)/ν11 at every position r, and additional DFTe energy components would modify this equilibrium. Predator–prey relationships are introduced by making species k a resource ({rho }_{k}=left]{n}_{k}right[), where (left]nright[) declares n a constant w.r.t. the functional differentiation of E, that is, the predator tends to align with the prey, not the prey with the predator. In view of the energy minimisation, the quadratic term in Eq. (7) entails that regions of low resource density ρk are less important than regions of high ρk. The different resources k have the same ability to limit the abundances, such that the limiting resource k = l at r has to come with the largest of weights wl(r), irrespective of the absolute amounts of resources at r. For example, the weights wk have to ensure that an essential but scarce mineral has (a priori) the same ability to limit the abundances as a resource like water, which might be abundant in absolute terms. To that end, we specify the weights$${w}_{k}({{{{{{{bf{r}}}}}}}})=frac{1}{{bar{rho }}_{k}^{2}}mathop{sum}limits_{s}eta ({nu }_{ks})exp left[sigma left(frac{{lambda }_{ks}}{{lambda }_{ls}}-1right)right],$$
(8)
which are inspired by the smooth minimum function, where σ  λls irrelevant at r. Using ({E}_{{{{{{{{rm{Res}}}}}}}}}), we show that an analytically solvable minimal example of two amensalistically interacting species already exhibits a plethora of resource-dependent equilibrium states (see Supplementary Notes and Supplementary Fig. 1).We specify the DFTe energy functional in Eq. (1) by summing Eqs. (2), (3), (6), and (7) and by (optionally) constraining the abundances to N via Lagrange multipliers μ:$$E[{{{{{{{bf{n}}}}}}}},, {{{{{{{boldsymbol{mu }}}}}}}}]({{{{{{{bf{N}}}}}}}}) equiv E[{{{{{{{bf{n}}}}}}}}]+{E}_{{{{{{{{boldsymbol{mu }}}}}}}}}[{{{{{{{bf{n}}}}}}}}]({{{{{{{bf{N}}}}}}}})\ equiv {E}_{{{{{{{{rm{dis}}}}}}}}}[{{{{{{{bf{n}}}}}}}}]+{E}_{{{{{{{{rm{env}}}}}}}}}[{{{{{{{bf{n}}}}}}}}]+{E}_{gamma }[{{{{{{{bf{n}}}}}}}}]+{E}_{{{{{{{{rm{Res}}}}}}}}}[{{{{{{{bf{n}}}}}}}}]+mathop{sum }limits_{s=1}^{S}{mu }_{s}left({N}_{s}-{int}_{A}({{{{{{{rm{d}}}}}}}}{{{{{{{bf{r}}}}}}}}),{n}_{s}right).$$
(9)
Uniform situations are characterised by spatially constant ingredients ns = Ns/A, ρk = Rk/A, coefficients τs, etc. for the DFTe energy, such that Eq. (9) reduces to a function E(N) with building blocks$${E}_{{{{{{{{rm{dis}}}}}}}}}longrightarrow frac{1}{2,A}mathop{sum }limits_{s=1}^{S}{tau }_{s},{N}_{s}^{2},$$
(10)
$${E}_{{{{{{{{rm{env}}}}}}}}}longrightarrow mathop{sum }limits_{s=1}^{S}{V}_{s}^{{{{{{{{rm{env}}}}}}}}},{N}_{s},$$
(11)
$${E}_{gamma }longrightarrow mathop{sum }limits_{{s,{s}^{{prime} }!=!1}atop {{s}^{{prime} }ne s}}^{S}frac{{N}_{s}^{{alpha }_{s}},{gamma }_{s{s}^{{prime} }},{N}_{{s}^{{prime} }}^{{beta }_{{s}^{{prime} }}}}{{A}^{{alpha }_{s}+{beta }_{{s}^{{prime} }}-1}},$$
(12)
$${E}_{{{{{{{{rm{Res}}}}}}}}}longrightarrow Amathop{sum }limits_{k=1}^{K}{{{{{{{{mathcal{L}}}}}}}}}_{k}left({{{{{{{bf{N}}}}}}}}/A,, {R}_{k}/Aright).$$
(13)
Ecosystem equilibria from the DFTe energy functionalThe general form of Eq. (9) gives rise to two types of minimisers (viz., equilibria): First, we term$${{{{{{{mathcal{H}}}}}}}}({{{{{{{bf{N}}}}}}}})equiv E[tilde{{{{{{{{bf{n}}}}}}}}}]equiv mathop{min }limits_{{{{{{{{bf{n}}}}}}}}}left{E[{{{{{{{bf{n}}}}}}}}],left|,{int}_{A}({{{{{{{rm{d}}}}}}}}{{{{{{{bf{r}}}}}}}}),{{{{{{{bf{n}}}}}}}}({{{{{{{bf{r}}}}}}}})={{{{{{{bf{N}}}}}}}},{{{{{{{rm{(fixed)}}}}}}}}right.right}$$
(14)
the ‘DFTe hypersurface’, with (tilde{{{{{{{{bf{n}}}}}}}}}) the energy-minimising spatial density profiles for given (fixed) N. Second, the ecosystem equilibrium is attained at the equilibrium abundances (hat{{{{{{{{bf{N}}}}}}}}}={int}_{A}({{{{{{{rm{d}}}}}}}}{{{{{{{bf{r}}}}}}}}),hat{{{{{{{{bf{n}}}}}}}}}({{{{{{{bf{r}}}}}}}})), which yield the global energy minimum$${{{{{{{mathcal{H}}}}}}}}(hat{{{{{{{{bf{N}}}}}}}}})=mathop{min }limits_{{{{{{{{bf{N}}}}}}}}},{{{{{{{mathcal{H}}}}}}}}({{{{{{{bf{N}}}}}}}}),$$
(15)
where the minimisation samples all admissible abundances, that is, ({{{{{{{bf{N}}}}}}}}in {left({{mathbb{R}}}_{0}^{+}right)}^{times S}) if no further constraints are imposed.The direct minimisation of E[n] is most practical for uniform systems, which only require us to minimise E(N) over an S-dimensional space of abundances. For the general nonuniform case, we adopt a two-step strategy that reflects Eqs. (14) and (15). First, we obtain the equilibrated density distributions on ({{{{{{{mathcal{H}}}}}}}}) for fixed N from the computational DPFT framework26,27,28,29,30,31. Second, a conjugate gradient descent searches ({{{{{{{mathcal{H}}}}}}}}({{{{{{{bf{N}}}}}}}})) for the global minimiser (hat{{{{{{{{bf{N}}}}}}}}}). Technically, we perform the computationally more efficient descent in μ-space. Local minima are frequently encountered, and we identify the best candidate for the global minimum from many individual runs that are initialised with random μ. Note that system realisations with energies close to the global minimum, especially local minima, are likely observable in reality, assuming that the system can equilibrate at all. There is always an equilibrium if the energy functional is bounded from below, together with the fact that the support (abundances/densities) of the energy functional is finite in any practical application. If some DFTe energy components are chosen (too) negative, the system can be unstable, in which case the energy functional has no minimum and is inappropriate for modelling the equilibrium in question. This means that another energy functional has to be considered, or, in the worst case, that DFTe is incapable of simulating this system. We also caution that no numerical optimisation algorithms for non-convex black-box functions can guarantee to find the global minimum, not even approximately. Without analytically available characteristics of the global minimum, all one may hope for are candidates of the minimiser, and those may not even be local minima—there is no way to be certain that an optimum proposed by a numerical optimisation algorithm is stable.Density-potential functional theory (DPFT) in Thomas–Fermi (TF) approximationDefining$${V}_{s}({{{{{{{bf{r}}}}}}}})={mu }_{s}-frac{delta {E}_{{{{{{{{rm{dis}}}}}}}}}[{{{{{{{bf{n}}}}}}}}]}{delta {n}_{s}({{{{{{{bf{r}}}}}}}})}$$
(16)
for all s, we obtain the reversible Legendre transform$${E}_{{{{{{{{rm{dis}}}}}}}}}^{{{{{{{{rm{L}}}}}}}}}[{{{{{{{bf{V}}}}}}}}-{{{{{{{boldsymbol{mu }}}}}}}}]={E}_{{{{{{{{rm{dis}}}}}}}}}[{{{{{{{bf{n}}}}}}}}]+mathop{sum }limits_{s=1}^{S}{int}_{A}({{{{{{{rm{d}}}}}}}}{{{{{{{bf{r}}}}}}}}),({V}_{s}-{mu }_{s}),{n}_{s}$$
(17)
of the dispersal energy and thereby supplement the total energy with the additional variables V:$$E[{{{{{{{bf{V}}}}}}}},, {{{{{{{bf{n}}}}}}}},, {{{{{{{boldsymbol{mu }}}}}}}}]({{{{{{{bf{N}}}}}}}})={E}_{{{{{{{{rm{dis}}}}}}}}}^{{{{{{{{rm{L}}}}}}}}}[{{{{{{{bf{V}}}}}}}}-{{{{{{{boldsymbol{mu }}}}}}}}]-{int}_{A}({{{{{{{rm{d}}}}}}}}{{{{{{{bf{r}}}}}}}}),{{{{{{{bf{n}}}}}}}}cdot ({{{{{{{bf{V}}}}}}}}-{{{{{{{{bf{V}}}}}}}}}^{{{{{{{{rm{env}}}}}}}}})+{E}_{{{{{{{{rm{int}}}}}}}}}[{{{{{{{bf{n}}}}}}}}]+{{{{{{{boldsymbol{mu }}}}}}}}cdot {{{{{{{bf{N}}}}}}}}.$$
(18)
This density-potential functional is equivalent to (but more flexible than) the density-only functional E[n,  μ](N). The minimisers of E[n] are thus among the stationary points of Eq. (18) and are obtained by solving$${n}_{s}[{V}_{s}-{mu }_{s}]({{{{{{{bf{r}}}}}}}})=frac{delta {E}_{{{{{{{{rm{dis}}}}}}}}}^{{{{{{{{rm{L}}}}}}}}}[{V}_{s}-{mu }_{s}]}{delta {V}_{s}({{{{{{{bf{r}}}}}}}})}$$
(19)
and$${V}_{s}[{{{{{{{bf{n}}}}}}}}]({{{{{{{bf{r}}}}}}}})={V}_{s}^{{{{{{{{rm{env}}}}}}}}}({{{{{{{bf{r}}}}}}}})+frac{delta {E}_{{{{{{{{rm{int}}}}}}}}}[{{{{{{{bf{n}}}}}}}}]}{delta {n}_{s}({{{{{{{bf{r}}}}}}}})}$$
(20)
self-consistently for all ns while enforcing ∫A(dr) ns(r) = Ns. Specifically, starting from V(0) = Venv, such that ({n}_{s}^{(0)}={n}_{s}[{V}_{s}^{(0)}-{mu }_{s}^{(0)}]), we iterate$${n}_{s}^{(i)}mathop{longrightarrow }limits^{{{{{{{{rm{equation}}}}}}}},(20)}{V}_{s}^{(i+1)}={V}_{s}[{{{{{{{{bf{n}}}}}}}}}^{(i)}]mathop{longrightarrow }limits^{{{{{{{{rm{equation}}}}}}}},(19)}{n}_{s}^{(i+1)}=(1-{theta }_{s}),{n}_{s}^{(i)}+{theta }_{s},{n}_{s}left[{V}_{s}^{(i+1)}-{mu }_{s}^{(i+1)}right]$$
(21)
until all ns are converged sufficiently. This self-consistent loop establishes a trade-off between dispersal energy and effective environment V by forcing an initial out-of-equilibrium density distribution to equilibrate at fixed N. We adjust ({mu }_{s}^{(i)}) in each iteration i such that ({n}_{s}^{(i)}) integrates to Ns. Small enough density admixtures, with 0  More

Martins, L. C. et al. Poluição atmosférica e atendimentos por pneumonia e gripe em São Paulo, Brasil. Revista de Saúde Pública 36, 88–94 (2002).Article 
PubMed 

Google Scholar 
Goudarzi, G. et al. Health risk assessment on human exposed to heavy metals in the ambient air PM10 in Ahvaz, Southwest Iran. Int. J. Biometeorol. 62, 1075–1083 (2018).Article 
ADS 
PubMed 

Google Scholar 
Makri, A. & Stilianakis, N. I. Vulnerability to air pollution health effects. Int. J. Hygiene Environ. Health 211, 326–336 (2008).Article 

Google Scholar 
Idani, E. et al. Characteristics, sources, and health risks of atmospheric PM10-bound heavy metals in a populated Middle Eastern City. Toxin Rev. 39, 266–274 (2020).Article 

Google Scholar 
Wang, J., Hu, Z., Chen, Y., Chen, Z. & Xu, S. Contamination characteristics and possible sources of PM10 and PM2.5 in different functional areas of Shanghai, China. Atmos. Environ. 68, 221–229 (2013).Article 
ADS 
CAS 

Google Scholar 
Guarnieri, M. & Balmes, J. R. Outdoor air pollution and asthma. Lancet 383, 1581–1592 (2014).Article 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Anderson, J. O., Thundiyil, J. G. & Stolbach, A. Clearing the air: A review of the effects of particulate matter air pollution on human health. J. Med. Toxicol. 8, 166–175 (2012).Article 
CAS 
PubMed 

Google Scholar 
Roy, D., Seo, Y.-C., Kim, S. & Oh, J. Human health risks assessment for airborne PM10-bound metals in Seoul, Korea. Environ. Sci. Pollut. Res. 26, 24247–24261 (2019).Article 
CAS 

Google Scholar 
Maesano, C. et al. Impacts on human mortality due to reductions in PM10 concentrations through different traffic scenarios in Paris, France. Sci. The Total. Environ. 698, 134257 (2020).Article 
ADS 
CAS 

Google Scholar 
Maleki, H., Sorooshian, A., Goudarzi, G., Nikfal, A. & Baneshi, M. M. Temporal profile of PM10 and associated health effects in one of the most polluted cities of the world (Ahvaz, Iran) between 2009 and 2014. Aeolian Res. 22, 135–140 (2016).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Medina, S., Le Tertre, A. & Saklad, M. The Apheis project: Air pollution and health—A European information system. Air Qual. Atmos. Heal. 2, 185–198 (2009).Article 

Google Scholar 
Medina, S., Plasencia, A., Ballester, F., Mücke, H. & Schwartz, J. Apheis: Public health impact of PM10 in 19 European cities. J. Epidemiol. Community Heal. 58, 831–836 (2004).Article 
CAS 

Google Scholar 
Pérez-Martínez, P. J., de Fátima Andrade, M. & de Miranda, R. M. Traffic-related air quality trends in São Paulo, Brazil. J. Geophys. Res. Atmos. 120, 6290–6304 (2015).Article 
ADS 

Google Scholar 
Sánchez-Ccoyllo, O. R. et al. Vehicular particulate matter emissions in road tunnels in Sao Paulo, Brazil. Environ. Monitoring Assess. 149, 241–249 (2009).Article 

Google Scholar 
Ribeiro, H. & de Assunção, J. V. Historical overview of air pollution in São Paulo Metropolitan Area, Brazil: Influence of mobile sources and related health effects. WIT Trans. Built Environ. 52,10 (2001).Bravo, M. A. & Bell, M. L. Spatial heterogeneity of PM10 and O3 in São Paulo, Brazil, and implications for human health studies. J. Air Waste Manag. Assoc. 61, 69–77 (2011).Article 
CAS 
PubMed 

Google Scholar 
De Freitas, E. D., Martins, L. D., da Silva Dias, P. L. & de Fátima Andrade, M. A simple photochemical module implemented in rams for tropospheric ozone concentration forecast in the metropolitan area of Sao Paulo, Brazil: Coupling and validation. Atmos. Environ. 39, 6352–6361 (2005).Article 
ADS 

Google Scholar 
Encalada-Malca, A. A., Cochachi-Bustamante, J. D., Rodrigues, P. C., Salas, R. & López-Gonzales, J. L. A spatio-temporal visualization approach of PM10 concentration data in Metropolitan Lima. Atmosphere 12, 609 (2021).Article 
ADS 

Google Scholar 
do Meio Ambiente, C. N. Institutes the national air quality control programee. Tech. Rep., Official Journal of the Federative Republic of Brazil (1989).do Meio Ambiente, C. N. Sets standards of primary and secondary air quality and even the criteria for acute episodes of air pollution. Tech. Rep., Official Journal of the Federative Republic of Brazil (1990).Artaxo, P. O estado da qualidade do ar no brasil. Work. Pap. WRI Brasil 32 (2021).Costa, A. F., Hoek, G., Brunekreef, B. & Ponce de Leon, A. C. Air pollution and deaths among elderly residents of Sao Paulo, Brazil: An analysis of mortality displacement. Environ. Health Perspectives 125, 349–354 (2017).Article 
CAS 

Google Scholar 
Bravo, M. A., Son, J., De Freitas, C. U., Gouveia, N. & Bell, M. L. Air pollution and mortality in São Paulo, Brazil: Effects of multiple pollutants and analysis of susceptible populations. J. Exposure Sci. Environ. Epidemiol. 26, 150–161 (2016).Article 
CAS 

Google Scholar 
Chiarelli, P. S. et al. The association between air pollution and blood pressure in traffic controllers in Santo André, São Paulo, Brazil. Environ. Res. 111, 650–655 (2011).Article 
CAS 

Google Scholar 
Ventura, L. M. B., de Oliveira Pinto, F., Soares, L. M., Luna, A. S. & Gioda, A. Forecast of daily PM2.5 concentrations applying artificial neural networks and holt-winters models. Air Qual. Atmos. Heal. 12, 317–325 (2019).Article 
CAS 

Google Scholar 
Leão, M. L. P., Zhang, L. & da Silva Júnior, F. M. R. Effect of particulate matter (PM2.5 and PM10) on health indicators: Climate change scenarios in a Brazilian Metropolis. Environ. Geochem. Heal. 44, 1–12 (2022).Habermann, M. & Gouveia, N. Application of land use regression to predict the concentration of inhalable particular matter in São Paulo City, Brazil. Engenharia Sanit. e Ambiental 17, 155–162 (2012).Article 

Google Scholar 
Braga, A. L. F., Pereira, L. A. A., Procópio, M., André, P. A. D. & Saldiva, P. H. D. N. Association between air pollution and respiratory and cardiovascular diseases in Itabira, Minas Gerais State. Brazil. Cadernos de Saúde Pública 23, S570–S578 (2007).Article 
PubMed 

Google Scholar 
Pinto, W. D. P., Reisen, V. A. & Monte, E. Z. Previsão da concentração de material particulado inalável, na região da grande vitória, ES, Brasil, utilizando o modelo sarimax. Engenharia Sanitária e Ambiental 23, 307–318 (2018).Article 

Google Scholar 
Schornobay-Lui, E. et al. Prediction of short and medium term PM10 concentration using artificial neural networks. Manag. Environ. Qual. An Int. J. 30, 414–436 (2018).Neto, P. S. D. M. et al. Neural-based ensembles for particulate matter forecasting. IEEE Access 9, 14470–14490 (2021).Article 

Google Scholar 
Albuquerque Filho, F. S. D., Madeiro, F., Fernandes, S. M., de Mattos Neto, P. S. & Ferreira, T. A. Time-series forecasting of pollutant concentration levels using particle swarm optimization and artificial neural networks. Química Nova 36, 783–789 (2013).Lei, T. M., Siu, S. W., Monjardino, J., Mendes, L. & Ferreira, F. Using machine learning methods to forecast air quality: A case study in Macao. Atmosphere 13, 1412 (2022).Article 
ADS 
CAS 

Google Scholar 
Yu, T. et al. Study on the regional prediction model of PM2.5 concentrations based on multi-source observations. Atmos. Pollut. Res. 13, 101363 (2022).Article 
CAS 

Google Scholar 
Li, J., Xu, G. & Cheng, X. Combining spatial pyramid pooling and long short-term memory network to predict PM2.5 concentration. Atmos. Pollut. Res. 13, 101309 (2022).Article 
CAS 

Google Scholar 
Cordova, C. H. et al. Air quality assessment and pollution forecasting using artificial neural networks in Metropolitan Lima-Peru. Sci. Rep. 11, 1–19 (2021).Article 

Google Scholar 
Plocoste, T., Calif, R. & Jacoby-Koaly, S. Temporal multiscaling characteristics of particulate matter PM10 and ground-level ozone O3 concentrations in caribbean region. Atmos. Environ. 169, 22–35 (2017).Article 
ADS 
CAS 

Google Scholar 
Calif, R. & Schmitt, F. G. Multiscaling and joint multiscaling description of the atmospheric wind speed and the aggregate power output from a wind farm. Nonlinear Process. Geophys. 21, 379–392 (2014).Article 
ADS 

Google Scholar 
Hyndman, R. J. & Khandakar, Y. Automatic time series forecasting: The forecast package for r. J. Stat. Softw. 27, 1–22 (2008).Article 

Google Scholar 
Harvey, A. C. Forecasting, structural time series models and the Kalman filter (Cambridge University Press, 1990).Book 
MATH 

Google Scholar 
Zhang, G. P. Time series forecasting using a hybrid arima and neural network model. Neurocomputing 50, 159–175 (2003).Article 
MATH 

Google Scholar 
Liao, T. W. Clustering of time series data—A survey. Pattern Recognit. 38, 1857–1874 (2005).Article 
ADS 
MATH 

Google Scholar 
Bell, M. L., Samet, J. M. & Dominici, F. Time-series studies of particulate matter. Annu. Rev. Public Heal. 25, 247–280 (2004).Article 

Google Scholar 
Hyndman, R. J. & Athanasopoulos, G. Forecasting: Principles and Practice (OTexts, 2018).
Google Scholar 
Box, G. E., Hillmer, S. C. & Tiao, G. C. Analysis and modeling of seasonal time series. in Seasonal analysis of economic time series, 309–344 (NBER, 1978).Sulandari, W., Suhartono, Subanar & Rodrigues, P. C. Exponential smoothing on modeling and forecasting multiple seasonal time series: An overview. Fluctuation Noise Lett. 20, 2130003 (2021).Article 
ADS 

Google Scholar 
Rodrigues, P. C., Awe, O. O., Pimentel, J. S. & Mahmoudvand, R. Modelling the behaviour of currency exchange rates with singular spectrum analysis and artificial neural networks. Stats 3, 137–157 (2020).Article 

Google Scholar 
Sako, K., Mpinda, B. N. & Rodrigues, P. C. Neural networks for financial time series forecasting. Entropy 24, 657 (2022).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Coelho, Leite et al. Statistical and artificial neural networks models for electricity consumption forecasting in the Brazilian industrial sector. Energies 15, 588 (2022).Article 

Google Scholar 
Sulandari, W., Subanar, S., Lee, M. H. & Rodrigues, P. C. Time series forecasting using singular spectrum analysis, fuzzy systems and neural networks. MethodsX 7, 101015 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Sulandari, W. et al. Indonesian electricity load forecasting using singular spectrum analysis, fuzzy systems and neural networks. Energy 190, 116408 (2020).Article 

Google Scholar 
Rodrigues, P. C. & Mahmoudvand, R. The benefits of multivariate singular spectrum analysis over the univariate version. J. Frankl. Inst. 355, 544–564 (2018).Article 
MathSciNet 
MATH 

Google Scholar  More

Different picocyanobacterial communities exhibit distinct gene repertoiresTo analyze the distribution of Prochlorococcus and Synechococcus reads along the Tara Oceans transect, metagenomic reads corresponding to the bacterial size fraction were recruited against 256 picocyanobacterial reference genomes, including SAGs and MAGs representative of uncultured lineages (e.g., Prochlorococcus HLIII-IV, Synechococcus EnvA or EnvB). This yielded a total of 1.07 billion recruited reads, of which 87.7% mapped onto Prochlorococcus genomes and 12.3% onto Synechococcus genomes, which were then functionally assigned by mapping them onto the manually curated Cyanorak v2.1 CLOG database [19]. In order to identify picocyanobacterial genes potentially involved in niche adaptation, we analyzed the distribution across the oceans of flexible (i.e. non-core) genes. Clustering of Tara Oceans stations according to the relative abundance of flexible genes resulted in three well-defined clusters for Prochlorococcus (Fig. 1A), which matched those obtained when stations were clustered according to the relative abundance of Prochlorococcus ESTUs, as assessed using the high-resolution marker gene petB, encoding cytochrome b6 (Fig. 1A; [24]). Only a few discrepancies can be observed between the two trees, including stations TARA-070 that displayed one of the most disparate ESTU compositions and TARA-094, dominated by the rare HLID ESTU (Fig. 1A). Similarly, for Synechococcus, most of the eight assemblages of stations discriminated based on the relative abundance of ESTUs (Fig. 1B) were also retrieved in the clustering based on flexible gene abundance, except for a few intra-assemblage switches between stations, notably those dominated by ESTU IIA (Fig. 1B). Despite these few variations, four major clusters can be clearly delineated in both Synechococcus trees, corresponding to four broadly defined ecological niches, namely (i) cold, nutrient-rich, pelagic or coastal environments (blue and light red in Fig. 1B), (ii) Fe-limited environments (purple and grey), (iii) temperate, P-depleted, Fe-replete areas (yellow) and (iv) warm, N-depleted, Fe-replete regions (dark red). This correspondence between taxonomic and functional information was also confirmed by the high congruence between distance matrices based on ESTU relative abundance and on CLOG relative abundance (p-value  0.01) are marked by a cross. Φsat: index of iron limitation derived from satellite data. PAR30: satellite-derived photosynthetically available radiation at the surface, averaged on 30 days. DCM: depth of the deep chlorophyll maximum.Full size imageIdentification of individual genes potentially involved in niche partitioningTo identify genes relevant for adaptation to a specific set of environmental conditions and enriched in specific ESTU assemblages, we selected the most representative genes from each module (Dataset 5; Figs. 3, S2). Most genes retrieved this way encode proteins of unknown or hypothetical function (85.7% of 7,485 genes). However, among the genes with a functional annotation (Dataset 6), a large fraction seems to have a function related to their realized environmental niche (Figs. 3, S2). For instance, many genes involved in the transport and assimilation of nitrite and nitrate (nirA, nirX, moaA-C, moaE, mobA, moeA, narB, M, nrtP; [6]) as well as cyanate, an organic form of nitrogen (cynA, B, D, S), are enriched in the Prochlorococcus blue module, which is correlated with the HLIIA-D ESTU and to low inorganic N, P, and silica levels and anti-correlated with Fe availability (Fig. 2A–C). This is consistent with previous studies showing that while only a few Prochlorococcus strains in culture possess the nirA gene and even less the narB gene, natural Prochlorococcus populations inhabiting N-poor areas do possess one or both of these genes [40,41,42]. Similarly, numerous genes amongst the most representative of Prochlorococcus brown, red and turquoise modules are related to adaptation of HLIIIA/IVA, HLIA and LLIA ESTUs to Fe-limited, cold P-limited, and cold, mixed waters, respectively (Fig. 3). Comparable results were obtained for Synechococcus, although the niche delineation was less clear than for Prochlorococcus since genes within each module exhibited lower correlations with the module eigenvalue (Fig. S2). These results therefore constitute a proof of concept that this network analysis was able to retrieve niche-related genes from metagenomics data.Fig. 3: Violin plots highlighting the most representative genes of each Prochlorococcus module.For each module, each gene is represented as a dot positioned according to its correlation with the eigengene for each module, the most representative genes being localized on top of each violin plot. Genes mentioned in the text and/or in Dataset 6 have been colored according to the color of the corresponding module, indicated by a colored bar above each module. The text above violin plots indicates the most significant environmental parameter(s) and/or ESTU(s) for each module, as derived from Fig. 2.Full size imageIdentification of eCAGs potentially involved in niche partitioningIn order to better understand the function of niche-related genes, notably of the numerous unknowns, we then integrated global distribution data with gene synteny in reference genomes using a network approach (Datasets 7, 8). This led us to identify clusters of adjacent genes in reference genomes, and thus potentially involved in the same metabolic pathway (Figs. 4, S3, S4; Dataset 6). These clusters were defined within each module and thus encompass genes with similar distribution and abundance in situ. Hereafter, these environmental clusters of adjacent genes will be called “eCAGs”.Fig. 4: Delineation of Prochlorococcus eCAGs, defined as a set of genes that are both adjacent in reference genomes and share a similar in situ distribution.Nodes correspond to individual genes with their gene name (or significant numbers of the CK number, e.g. 1234 for CK_00001234) and are colored according to their WGCNA module. A link between two nodes indicates that these two genes are less than five genes apart in at least one genome. The bottom insert shows the most significant environmental parameter(s) and/or ESTU(s) for each module, as derived from Fig. 2.Full size imageeCAGs related to nitrogen metabolismThe well-known nitrate/nitrite gene cluster involved in uptake and assimilation of inorganic forms of N (see above), which is present in most Synechococcus genomes (Dataset 6), was expectedly not restricted to a particular niche in natural Synechococcus populations, as shown by its quasi-absence from WGCNA modules. In Prochlorococcus, this cluster is separated into two eCAGs enriched in low-N areas (Fig. S5A, B), most genes being included in Pro-eCAG_002, present in only 13 out of 118 Prochlorococcus genomes, while nirA and nirX form an independent eCAG (Pro-eCAG_001) due to their presence in many more genomes. The quasi-core ureA-G/urtB-E genomic region was also found to form a Prochlorococcus eCAG (Pro-eCAG_003) that was impoverished in low-Fe compared to other regions (Fig. S5C, D), in agreement with its presence in only two out of six HLIII/IV genomes. We also uncovered several other Prochlorococcus and Synechococcus eCAGs that seem to be involved in the transport and/or assimilation of more unusual and/or complex forms of nitrogen, which might either be degraded into elementary N molecules or possibly directly used by cells for e.g. the biosynthesis of proteins or DNA. Indeed, we detected in both genera an eCAG (Pro-eCAG_004 and Syn-eCAG_001; Fig. S6A, B; Dataset 6) that encompasses speB2, an ortholog of Synechocystis PCC 6803 sll1077, previously annotated as encoding an agmatinase [29, 43] and which was recently characterized as a guanidinase that degrades guanidine rather than agmatine to urea and ammonium [44]. E. coli produces guanidine under nutrient-poor conditions, suggesting that guanidine metabolism is biologically significant and potentially prevalent in natural environments [44, 45]. Furthermore, the ykkC riboswitch candidate, which was shown to specifically sense guanidine and to control the expression of a variety of genes involved in either guanidine metabolism or nitrate, sulfate, or bicarbonate transport, is located immediately upstream of this eCAG in Synechococcus reference genomes, all genes of this cluster being predicted by RegPrecise 3.0 to be regulated by this riboswitch (Fig. S6C; [45, 46]). The presence of hypA and B homologs within this eCAG furthermore suggests that, in the presence of guanidine, these homologs could be involved in the insertion of Ni2+, or another metal cofactor, in the active site of guanidinase. The next three genes of this eCAG, which encode an ABC transporter similar to the TauABC taurine transporter in E. coli (Fig. S6C), could be involved in guanidine transport in low-N areas. Of note, the presence in most Synechococcus/Cyanobium genomes possessing this eCAG of a gene encoding a putative Rieske Fe-sulfur protein (CK_00002251) downstream of this gene cluster, seems to constitute a specificity compared to the homologous gene cluster in Synechocystis sp. PCC 6803. The presence of this Fe-S protein suggests that Fe is used as a cofactor in this system and might explain why this gene cluster is absent from picocyanobacteria thriving in low-Fe areas, while it is present in a large proportion of the population in most other oceanic areas (Fig. S6A, B).Another example of the use of organic N forms concerns compounds containing a cyano radical (C ≡ N). The cyanate transporter genes (cynABD) were indeed found in a Prochlorococcus eCAG (Pro-eCAG_005, also including the conserved hypothetical gene CK_00055128; Fig. S7A, B). While only a small proportion of the Prochlorococcus community possesses this eCAG in warm, Fe-replete waters, it is absent from other oceanic areas in accordance with its low frequency in Prochlorococcus genomes (present in only two HLI and five HLII genomes). In Synechococcus these genes were not included in a module, and thus are not in an eCAG (Dataset 6; Fig. S7C), but seem widely distributed despite their presence in only a few Synechococcus genomes (mostly in clade III strains; [6, 47, 48]). Interestingly, we also uncovered a 7-gene eCAG (Pro-eCAG_006 and Syn-eCAG_002), encompassing a putative nitrilase gene (nitC), which also suggests that most Synechococcus cells and a more variable fraction of the Prochlorococcus population could use nitriles or cyanides in warm, Fe-replete waters and more particularly in low-N areas such as the Indian Ocean (Fig. 5A, B). The whole operon (nitHBCDEFG; Fig. 5C), called Nit1C, was shown to be upregulated in the presence of cyanide and to trigger an increase in the rate of ammonia accumulation in the heterotrophic bacterium Pseudomonas fluorescens [49], suggesting that like cyanate, cyanide could constitute an alternative nitrogen source in marine picocyanobacteria as well. However, given the potential toxicity of these C ≡ N-containing compounds [50], we cannot exclude that these eCAGs could also be devoted to cell detoxification [45, 47]. Such an example of detoxification has been described for arsenate and chromate that, as analogs of phosphate and sulfate respectively, are toxic to marine phytoplankton and must be actively exported out of the cells [51, 52].Fig. 5: Global distribution map of the eCAG involved in nitrile or cyanide transport and assimilation.A Prochlorococcus Pro-eCAG_006. B Synechococcus Syn-eCAG_002. C The genomic region in Prochlorococcus marinus MIT9301. The size of the circle is proportional to relative abundance of each genus as estimated based on the single-copy core gene petB and this gene was also used to estimate the relative abundance of other genes in the population. Black dots represent Tara Oceans stations for which Prochlorococcus or Synechococcus read abundance was too low to reach the threshold limit.Full size imageWe detected the presence of an eCAG encompassing asnB, pyrB2, and pydC (Pro-eCAG_007, Syn-eCAG_003, Fig. S8), which could contribute to an alternative pyrimidine biosynthesis pathway and thus provide another way for cells to recycle complex nitrogen forms. While this eCAG is found in only one fifth of HLII genomes and in quite specific locations for Prochlorococcus, notably in the Red Sea, it is found in most Synechococcus cells in warm, Fe-replete, N and P-depleted niches, consistent with its phyletic pattern showing its absence only from most clade I, IV, CRD1, and EnvB genomes (Fig. S8; Dataset 6). More generally, most N-uptake and assimilation genes in both genera were specifically absent from Fe-depleted areas, including the nirA/narB eCAG for Prochlorococcus, as mentioned by Kent et al. [36] as well as guanidinase and nitrilase eCAGs. In contrast, picocyanobacterial populations present in low-Fe areas possess, in addition to the core ammonium transporter amt1, a second transporter amt2, also present in cold areas for Synechococcus (Fig. S9). Additionally, Prochlorococcus populations thriving in HNLC areas also possess two amino acid-related eCAGs that are present in most Synechococcus genomes, the first one involved in polar amino acid N-II transport (Pro-eCAG_008; natF-G-H-bgtA; [53]; Fig. S10A, B) and the second one (leuDH-soxA-CK_00001744, Pro-eCAG_009, Fig. S10C, D) that notably encompasses a leucine dehydrogenase, able to produce ammonium from branched-chain amino acids. This highlights the profound difference in N acquisition mechanisms between HNLC regions and Fe-replete, N-deprived areas: the primary nitrogen sources for picocyanobacterial populations dwelling in HNLC areas seem to be ammonium and amino acids, while N acquisition mechanisms are more diverse in N-limited, Fe-replete regions.eCAGs related to phosphorus metabolismAdaptation to P depletion has been well documented in marine picocyanobacteria showing that while in P-replete waters Prochlorococcus and Synechococcus essentially rely on inorganic phosphate acquired by core transporters (PstSABC), strains isolated from low-P regions and natural populations thriving in these areas additionally contain a number of accessory genes related to P metabolism, located in specific genomic islands [6, 14, 30,31,32, 54]. Here, we indeed found in Prochlorococcus an eCAG containing the phoBR operon (Pro-eCAG_010) that encodes a two-component system response regulator, as well as an eCAG including the alkaline phosphatase phoA (Pro-eCAG_011), both present in virtually the whole Prochlorococcus population from the Mediterranean Sea, the Gulf of Mexico and the Western North Atlantic Ocean, which are known to be P-limited [30, 55] (Fig. S11A, B). By comparison, in Synechococcus, we only identified the phoBR eCAG (Syn-eCAG_005, Fig. S11C) that is systematically present in warm waters whatever the limiting nutrient, in agreement with its phyletic pattern in reference genomes showing its specific absence from cold thermotypes (clades I and IV, Dataset 6). Furthermore, although our analysis did not retrieve them within eCAGs due to the variability of gene content and synteny in this genomic region, even within each genus, several other P-related genes were enriched in low-P areas but partially differed between Prochlorococcus and Synechococcus (Figs. 3, S2, S11; Dataset 6). While the genes putatively encoding a chromate transporter (ChrA) and an arsenate efflux pump ArsB were present in both genera in different proportions, a putative transcriptional phosphate regulator related to PtrA (CK_00056804; [56]) was specific to Prochlorococcus. Synechococcus in contrast harbors a large variety of alkaline phosphatases (PhoX, CK_00005263 and CK_00040198) as well as the phosphate transporter SphX (Fig. S11).Phosphonates, i.e. reduced organophosphorus compounds containing C–P bonds that represent up to 25% of the high-molecular-weight dissolved organic P pool in the open ocean, constitute an alternative P form for marine picocyanobacteria [57]. We indeed identified, in addition to the core phosphonate ABC transporter (phnD1-C1-E1), a second previously unreported putative phosphonate transporter phnC2-D2-E2-E3 (Pro-eCAG_012; Fig. 6A). Most of the Prochlorococcus population in strongly P-limited areas of the ocean harbored these genes, while they were absent from other areas, consistent with their presence in only a few Prochlorococcus and no Synechococcus genomes. Furthermore, as previously described [58,59,60], we found a Prochlorococcus eCAG encompassing the phnYZ operon involved in C-P bond cleavage, the putative phosphite dehydrogenase ptxD, and the phosphite and methylphosphonate transporter ptxABC (Pro-eCAG_0013, Dataset 6; Fig. 6B, [60,61,62]). Compared to these previous studies that mainly reported the presence of these genes in Prochlorococcus cells from the North Atlantic Ocean, here we show that they actually occur in a much larger geographic area, including the Mediterranean Sea, the Gulf of Mexico, and the ALOHA station (TARA_132) in the North Pacific, even though they were present in a fairly low fraction of Prochlorococcus cells. These genes occurred in an even larger proportion of the Synechococcus population, although not found in an eCAG for this genus (Fig. S12; Dataset 6). Synechococcus cells from the Mediterranean Sea, a P-limited area dominated by clade III [24], seem to lack phnYZ, in agreement with the phyletic pattern of these genes in reference genomes, showing the absence of this two-gene operon in the sole clade III strain that possesses the ptxABDC gene cluster. In contrast, the presence of the complete gene set (ptxABDC-phnYZ) in the North Atlantic, at the entrance of the Mediterranean Sea, and in several clade II reference genomes rather suggests that it is primarily attributable to this clade. Altogether, our data indicate that part of the natural populations of both Prochlorococcus and Synechococcus would be able to assimilate phosphonate and phosphite as alternative P-sources in low-P areas using the ptxABDC-phnYZ operon. Yet, the fact that no picocyanobacterial genome except P. marinus RS01 (Fig. 6C) possesses both phnC2-D2-E2-E3 and phnYZ, suggests that the phosphonate taken up by the phnC2-D2-E2-E3 transporter could be incorporated into cell surface phosphonoglycoproteins that may act to mitigate cell mortality by grazing and viral lysis, as recently suggested [63].Fig. 6: Global distribution map of eCAGs putatively involved in phosphonate and phosphite transport and assimilation.A Prochlorococcus Pro-eCAG_012 putatively involved in phosphonate transport. B Prochlorococcus Pro-eCAG_013, involved in phosphonate/phosphite uptake and assimilation and phosphonate C-P bond cleavage. C The genomic region encompassing both phnC2-D2-E2-E3 and ptxABDC-phnYZ specific to P. marinus RS01. The size of the circle is proportional to relative abundance of Prochlorococcus as estimated based on the single-copy core gene petB and this gene was also used to estimate the relative abundance of other genes in the population. Black dots represent Tara Oceans stations for which Prochlorococcus read abundance was too low to reach the threshold limit.Full size imageeCAGs related to iron metabolismAs for macronutrients, it has been hypothesized that the survival of marine picocyanobacteria in low-Fe regions was made possible through several strategies, including the loss of genes encoding proteins that contain Fe as a cofactor, the replacement of Fe by another metal cofactor, and the acquisition of genes involved in Fe uptake and storage [14, 15, 36, 39, 64]. Accordingly, several eCAGs encompassing genes encoding proteins interacting with Fe were found in modules anti-correlated to HNLC regions in both genera. These include three subunits of the (photo)respiratory complex succinate dehydrogenase (SdhABC, Pro-eCAG_014, Syn-eCAG_006, Fig. S13; [65]) and Fe-containing proteins encoded in most abovementioned eCAGs involved in N or P metabolism, such as the guanidinase (Fig. S6), the NitC1 (Fig. 5), the pyrB2 (Fig. S8), the phosphonate (Fig. 6, S12), and the urea and inorganic nitrogen eCAGs (Fig. S5). Most Synechococcus cells thriving in Fe-replete areas also possess the sodT/sodX eCAG (Syn-eCAG_007, Fig. S14A, B) involved in nickel transport and maturation of the Ni-superoxide dismutase (SodN), these three genes being in contrast core in Prochlorococcus. Additionally, Synechococcus from Fe-replete areas, notably from the Mediterranean Sea and the Indian Ocean, specifically possess two eCAGs (Syn-eCAG_008 and 009; Fig. S14C, D), involved in the biosynthesis of a polysaccharide capsule that appear to be most similar to the E. coli groups 2 and 3 kps loci [66]. These extracellular structures, known to provide protection against biotic or abiotic stress, were recently shown in Klebsiella to provide a clear fitness advantage in nutrient-poor conditions since they were associated with increased growth rates and population yields [67]. However, while these authors suggested that capsules may play a role in Fe uptake, the significant reduction in the relative abundance of kps genes in low-Fe compared to Fe-replete areas (t-test p-value  More

Genome reduction is associated with a termite comb-associated lifestyleFor our studies, we collected fungus comb samples originating from mounds of Macrotermes natalensis, Odontotermes spp., and Microtermes spp. termites and were able to obtain seven viable Pseudoxylaria cultures (X802 [Microtermes sp.], Mn132, Mn153, X187, X3-2 [Macrotermes natalensis], and X167, X170LB [Odontotermes spp.], Table S1-S3).To test if a fungus comb-associated lifestyle of Pseudoxylaria was reflected in differences at the genome level, we sequenced the genomes of all seven isolates using a combination of paired-end shotgun sequencing (BGISEQ-500, BGI) and long-read sequencing (PacBio sequel, BGI or Oxford Nanopore Technologies, Oxford, UK). In addition, we sequenced the transcriptomes (BGISEQ, BGI) of two isolates (X802, X170LB). Eleven publicly available genomes of free-living Xylaria (Fig. 2A, B) were used as reference genomes (Table S4). Hybrid draft genomes were comprised on average of 33–742 scaffolds with total haploid assembly lengths of 33.2–40.4 Mb, and a high BUSCO completeness of genomes ( > 95 %) with a total number of predicted proteins ranging from 8.8 to 12.1 × 103. The GC content was comparable to reference genomes with 49.7–51.6%. To verify the phylogenetic placement of the isolates, different genetic loci encoding conserved protein sequences (α-actin (ACT), second largest subunit of RNA polymerase (RPB2), β-tubulin (TUB) and the internal transcribed spacer (ITS) were used as genetic markers [7, 13].Fig. 2: Geographic and comparative phylogenomic analysis of termite-associated Pseudoxylaria isolates (strains 1-7) and free-living Xylaria (strains 8–18).A Geographic origins of genome-sequenced free-living Xylaria and termite-associated Pseudoxylaria isolates, B phylogenomic placement based on single-copy ortholog protein sequences, and C comparison of genome assembly length, and numbers of predicted proteins per genome.Full size imagePhylogenies were reconstructed from ITS sequences and three aligned sequence datasets (RPB2, TUB, ACT) using reference sequences of twelve different taxa (Table S4–S7). Consistent with previous findings, all isolates grouped within the monophyletic termite-associated Pseudoxylaria group [9,10,11,12,13], which diverged from the free-living members of the genus Xylaria (Fig. 2B, Figure S1–S4).As our seven isolates covered a larger portion of the previously reported phylogenetic diversity of the termite-associated subgenus, we elaborated on genomic characteristics of our isolates to uncover features of the termite-associated ecology of Pseudoxylaria. Indeed, comparative genome analysis of the South African Pseudoxylaria isolates with publicly available genomes of free-living Xylaria species of similar genome quality revealed significantly reduced genome assembly lengths in Pseudoxylaria with reduced numbers of predicted genes per genome (Table S4). Comparison of the annotated mitochondrial (mt) genomes (Figure S5, Table S8) also indicated that all seven mt genomes were shorter in length (assembly lengths: 18.5–63.8 kbp) compared to the, albeit few, publicly available mitochondrial genomes of free-living species (48.9–258.9 kbp). The reduction in mitochondrial genome size also corresponded to a significantly reduced mean number of annotated genes (7.6) and tRNAs (14.3) in Pseudoxylaria spp. compared to on average 30.0 (annotated genes) and 25.8 (tRNAs) found in free-living species.Analysis of the abundance and composition of transposable elements (TEs), which account for up to 30–35% of the genomes of (endo)parasitic fungi due to the expansion of certain gene families [20, 21], showed that the mean total numbers of TEs across Pseudoxylaria spp. genomes were comparable (1530), but the numbers were reduced compared to free-living Xylaria species (3690) (Table S9). We also identified high variation in the TE composition across genomes (1.5–9.9 %), comparable to what was observed in free-living Xylaria spp. (1.3–8.1 %), with reductions in long terminal repeat retrotransposons (LTRs: Copia and unknown LTRs) in two inverted tandem repeat DNA transposons (TIRs; CACTA, Mutator and hAT). As Pseudoxylaria spp. contained increased numbers of non-ITR transposons of the helitron class and LTRs of the Gypsy class compared to Xylaria strains, we concluded that Pseudoxylaria exhibits no typical features of an (endo)parasitic lifestyle, but that the overall composition and the reduced numbers of TEs could serve as a fingerprint to distinguish the genetically divergent Pseudoxylaria taxa.Repertoire of carbohydrate-active enzymes indicates specialized substrate useAs the fungus comb is mostly composed of partially-digested plant material interspersed with fungal mycelium of the termite mutualist [3], we anticipated that Pseudoxylaria should exhibit features of a substrate specialist similar to the fungal mutualist Termitomyces, which should be reflected in a Carbohydrate-Active enzyme (CAZyme) repertoire distinguishable from  free-living saprophytic Xylaria species [22,23,24]. In particular, numbers and composition of redox-active enzymes (e.g., benzoquinone reductase (EC 1.6.5.6/EC 1.6.5.7), catalase (EC 1.11.1.6), glutathione reductase (EC 1.11.1.9), hydroxy acid oxidase (EC 1.1.3.15), laccase (EC 1.10.3.2), manganese peroxidase (EC 1.11.1.13), peroxiredoxin (EC 1.11.1.15), superoxide dismutase (EC 1.15.1.1), dye-decolorization or unspecific peroxygenase (EC 1.11.2.1), Table S10), which catalyze the degradation of lignin-rich biomass, were expected to differ between free-living strains and substrate specialists [22].Identification of CAZymes using Peptide Pattern Recognition (PPR) revealed that Pseudoxylaria genomes encoded on average a reduced number of CAZymes (mean 264) compared to the free-living taxa in the family Xylaria (mean 367 CAZymes, pANOVA; F = 41.4, p = 3.5 × 10–8, pairwise p = 1.69 × 10–7) (Fig. 3A, B, Figure S6), but similar numbers to those identified in Termitomyces (mean 265, pairwise p = 0.949).Fig. 3: Comparison of carbohydrate-active enzymes (CAZymes) in Xylaria, Pseudoxylaria and the fungal mutualist Termitomyces.A Predicted CAZymes, B Principal Coordinates Analysis (PCoA) of predicted CAZyme families, and C heatmap of representatives CAZyme families in the predicted proteomes of free-living Xylaria, Termitomyces and Pseudoxylaria species.Full size imageOverall, significant differences in the composition of CAZymes were observed [8], most notably in the reduction of auxiliary activity enzymes (AA), carbohydrate esterases (CE), glycosyl hydrolases (GH), and polysaccharide lyases (PL). The most significant reduction was observed in the AA3 family (Fig. 3C), which typically displays a high multigenicity in wood-degrading fungi as many  enzymes of this family catalyze the oxidation of alcohols or carbohydrates with the concomitant formation of hydrogen peroxide or hydroquinones thereby supporting lignocellulose degradation by other AA-enzymes, such as peroxidases (AA2). Similarly, although to a lesser extent, reduced numbers within the related AA1 family were detected, which included oxidizing enzymes like laccases, ferroxidases, and laccase-like multicopper oxidases. Along these lines, glycosyl hydrolases of the GH3 and GH5 family, including enzymes responsible for degradation of cellulose-containing biomass and xylose, were less abundant. We also noted that all Pseudoxylaria lacked homologs of the unspecific peroxygenases (UPO; EC 1.11.2.1), while almost all free-living Xylaria spp. and the fungal symbiont Termitomyces harbored at least one or two copies of similar gene sequences.
Pseudoxylaria shows reduced biosynthetic capacity for secondary metabolite productionA healthy termite colony is engulfed in several layers of social immunity [5, 6], which pose a constant selection pressure on associated and potentially antagonistic microbes. As Pseudoxylaria evolved measures to remain inconspicuously present within the comb environment, we hypothesized that one of the possible adaptations to evade hygiene measures of termites could be reflected in a reduced biosynthetic capability to produce antibiotic or volatile natural products, which often serve as infochemicals triggering defense mechanisms [25,26,27], or as alarm pheromones [4, 28].The biosynthesis of secondary metabolites is encoded in so called Biosynthetic Gene Cluster (BGC) regions. We explored the abundance and diversity of encoded BGCs using FungiSMASH 6.0.0 and manually cross-checked the obtained data set by BLAST to account for possible biases due to varying genome qualities across strains of both groups [29]. Overall, the herein investigated Xylaria genomes harbored on average 90 BGCs per genome, while Pseudoxylaria encoded on average 45 BGCs (Fig. 4, Figure S7). Fig. 4: Similarity network analysis of biosynthetic gene clusters.Comparative analysis of termite associated-associated Pseudoxylaria isolates (strains 1–7, red circles) and free-living Xylaria (strains 8–18, green circles) with BiG-SCAPE 1.0 annotations (blue hexagon) ACR ACR toxin, Alt alternariol, Bio biotin, Chr chromene, Cyt cytochalasins, Cur curvupalide, Dep depiudecin, Fus fusarin, Gri griseofulvin, Mon monascorubin, MSA 6-methylsalicylic acid, Pho phomasetin, Sol solanapyrone, Swa swasionine, Xen xenolozoyenone, Xsp xylasporins, Xyl xylacremolide. Singletons are not shown.Full size imageThe nature and relatedness of the BGCs were analyzed by creating a curated similarity network analysis using BiG-SCAPE 1.0 [30]. Overall, 28 orthologous BGCs were shared across all genomes, including the biosynthesis of polyketides like 6-methylsalicylic acid (MSA), chromenes (Chr) and polyketide-non-ribosomal peptide (PKS-NRPS) hybrids like the cytochalasins (Cyt) [31]. Furthermore, five BGC networks, which were shared by Pseudoxylaria and Xylaria, contained genes encoding natural product modifying dimethylallyltryptophan synthases (DMATS). In contrast, and despite the significant reduction in the biosynthetic capacity within Pseudoxylaria genomes [29], about 29 BGC networks were unique to Pseudoxylaria and thus could possibly relate to the comb-associated lifestyle (Figure S8 and S9). Notably, Pseudoxylaria genomes lacked genes encoding ribosomally synthesized and posttranslationally modified peptides (RiPPs) or halogenases. In comparision, free-living Xylaria spp. harbored at least one sequence encoding a RiPP, and up to two orthologous sequences encoding putative halogenases. In contrast, a reduced average number of terpene synthases (TPS) in Pseudoxylaria (9 TPS) compared to free-living Xylaria (18 TPS) was detected, which included three BGCs encoding TPSs that were unique to Pseudoxylaria.  In comparison, genomes of the fungal mutualist Termitomyces were reported to encode for about 20-25 terpene cyclases, but haboured only about two loci containing genes for a PKS and NRPS each [24].Manual BLAST searches were conducted to identify BGCs that could be putatively assigned to previously isolated metabolites from Pseudoxylaria (vide infra Fig. 7, Figure S8) [32, 33]. Using e.g., the known NRPS-PKS-hybrid cluster sequence ccs (Aspergillus clavatus) of cytochalasins as query, an orthologous BGC, here named cytA, was identified in the cytochalasin-producing strain X802 [34]. Although the putative PKS-NRPS hybrid and CcsA shared 60 % identical amino acids (aa), the sequences of the accessory enzymes were less related to CcsC-G (45–47% identical aa) and the BGC in X802 lacked a gene of a homologue to ccsB. Similarly, five free-living Xylaria species carried orthologous gene loci (Xylaria sp. BCC 1067, Xylaria sp. MSU_SB201401, X. flabelliformis G536, X. grammica EL000614, and X. multiplex DSM 110363) supporting previous isolation reports of cytochalasins with varying structural features. Furthermore, three Pseudoxylaria strains (X187, and closely related Mn153, and Mn132) were found to share a highly similar PKS-NRPS hybrid BGC (99–100 % identical aa, named xya), which likely encodes for the enzymatic production of previously identified xylacremolides [32]. Four Pseudoxylaria strains (X802, Mn132, Mn153, and X187) also shared a BGC (50–98 % amino acid identity) resembling the fog BGC (Aspergillus ruber) [35, 36], which putatively encodes the biosynthetic machinery to produce xylasporin/cytosporin-like metabolites. In this homology search, we also uncovered that fog-like BGC arrangements are likely more common than previously anticipated, as clusters with similar arrangements and identity were also found in genomes of Rosellinia necatrix, Pseudomasariella vexata, Stachybotrys chartarum, and Hyaloscypha bicolor (Fig. 4, Figure S8).A detailed analysis of the fog-like cluster arrangements within Pseudoxylaria genomes revealed – similar to homologs of the ccs cluster – variation in the abundance and arrangement of several accessory genes coding for a cupin protein (pxF), a short chain oxidoreductase (pxB; SDR), and an additional SnoaL-like polyketide cyclase (pxP), which could account for the production of strain-specific structural congeners (vide infra, Fig. 7).Change of nutrient sources causes dedicated transcriptomic changes in Pseudoxylaria
To further solidify our in silico indications of substrate specialization with comb material as preferred substrate and fungus garden as environment, we analyzed Pseudoxylaria growth on different media (PDA, and reduced medium 1/3-PDA) including comb-like agar matrices (wood-rice medium (WRM), agar-agar or 1/3-PDA medium containing lyophilized (dead) Termitomyces sp. T112 biomass (T112, respectively T112-PDA), PDB covering glass-based surface-structuring elements (GB), Table S11–S14).Cultivation of Pseudoxylaria on agar-agar containing lyophilized biomass of Termitomyces (T112) as the sole nutrient source allowed Pseudoxylaria to sustain growth, although to a reduced extent compared to growth on nutrient-rich PDA medium (Table S3). Wood-rice medium (WRM) induced comparable growth rates as observed on PDA and also the appearance of phenotypic stromata.To investigate the influence of these growth conditions on the transcriptomic level, we harvested RNA from vegetative mycelium after growth on comb-like media (WRM, T112, T112-PDA, and GB), PDA, and reduced medium 1/3-PDA (Fig. 5A). The most significant transcript changes (normalized to data obtained from growth on PDA) were observed for genes coding for specific CAZymes including several redox active enzymes (Fig. 5B). The 30 most variable transcripts coded for specific glycoside hydrolases (GH), lytic polysaccharide monooxygenases (AA), ligninolytic enzymes, and a glycoside transferase (GT). Similarly, chitinases (CHT2; CHT4; CHI2; CHI4) were upregulated (up to 243-fold on T112) under almost all conditions compared to PDA, but some of these specific transcript changes were exclusive to growth on Termitomyces biomass or artificial comb material (WRM) suggesting the ability to regulate and increase chitin metabolism if necessary [37].Fig. 5: Transcriptomic analysis of Pseudoxylaria sp. X802 in dependence of growth conditions.A Representative pictures of Pseudoxylaria sp. X802 growing on PDA, PDB on glass beads (GB), wood-rice medium (WRM), and agar-agar medium containing lyophilized Termitomyces sp. T112 biomass (T112). B Heatmap of the most variable transcripts coding for CAZymes (red), redox enzymes (orange), secondary metabolite-related core genes (green), and more specifically on key genes within the boundaries of cytochalasin (turquoise) and xylasporin/cytosporin BGCs (blue). RNA was obtained from vegetative mycelium after growth on PDA, reduced medium (1/3-PDA), PDB on glass beads (GB), wood-rice medium (WRM), 1/3-PDA-medium enriched with Termitomyces sp. T112 biomass (T112-PDA) and agar-agar medium containing lyophilized Termitomyces biomass (T112). Transcript counts are shown as log10 transformed transcripts per million (top; TPM). Significance of the changes in transcript counts are compared to control (X802 grown on PDA) and depicted in log-10 transformed p values.Full size imageWhen X802 was grown on T112 (agar matrix containing lyophilized Termitomyces sp. T112 biomass), we observed a >400-fold increase in the expression of transcripts encoding glycoside hydrolases in the GH43 family, GH7 (~140-fold), GH3, and GH64 (5–12-fold). Similarly, transcripts for a putative mannosyl-oligosaccharide-α-1,2-mannosidase (MNS1B; 8.2-fold), chitinase CHT4 (2.9-fold), β-glucosidase BGL4 (5.7-fold), and copper-dependent lytic polysaccharide monooxygenase AA11 (1.6-fold) were significantly upregulated. Growth on WRM (wood-rice medium) or T112 (Termitomyces sp. T112 biomass) also caused a significant upregulation of genes coding for glycoside transferase GT2, glycoside hydrolases GH15, GH3, and aldehyde oxidase AOX1, which indicated the ability to expand the degradation portfolio if necessary. Along these lines, specific transcript levels were reduced when X802 was grown on T112, in particular class II lignin-modifying peroxidases (AA2), carbohydrate-binding module family 21 (CBM21), multicopper oxidases (AA1), secreted β-glucosidases (SUN4), and glycoside hydrolases GH16, and GH128.When the fungus was challenged with lignocellulose-rich WRM medium, higher transcript levels putatively assigned to glutathione peroxidase (GXP2), superoxide dismutase (SOD2), and laccases (LCC5) were observed, which indicated that despite the reduced wood-degrading capacity, Pseudoxylaria activates available enzymatic mechanisms to degrade the provided material and respond to the resulting oxidative stress. Cultivation on GB (glass-based surfaces covered in liquid PD broth) influenced the expression of certain genes coding for glycoside hydrolases (GH64, GH76, GH72, GH128, BGL4) and lytic polysaccharide monooxygenases (AA1, AA2, AA11), presumably enabling the fungus to utilize soluble carbohydrates.To test the hypothesis that the presence of Termitomyces biomass stimulates secondary metabolite production in Pseudoxylaria to eventually displace the mutualist, we also analyzed changes in the transcript levels of core BGC genes that encode the production of bioactive secondary metabolites. Overall, only slight transcript variations were detectable within the  most variable expressed genes. (Fig. 5B). Cultivation on GB, WRM, and T112 media caused lower transcript levels of genes coding for terpene synthase TC1, polyketide synthases (PKS7, PKS8), and the NRPS-like1, while an upregulation of NRPS-like2 on WRM (2.5-fold), and of PKS7 (1.7-fold) on reduced 1/3-PDA medium was observed.Transcript levels of core genes within BGCs assigned to cytochalasines (cyt) or xylasporins/cytosporins (px), e.g., remained nearly constant, while minor transcript level variations of neighboring genes and reduced transcript levels for pxI (flavin-dependent monooxygenase), pxH (ABBA-type prenyltransferase), pxF (cupin fold oxidoreductase), and pxJ (short-chain dehydrogenase) were detectable. Hence, it was concluded that the presence of Termitomyces biomass only weakly triggers secondary metabolite production in general, but varying transcript levels coding for decorating enzymes could cause substantial structural alterations within the produced natural product composition. It was also notable that transcript levels of the terpene synthase TC1 were downregulated, which could cause a reduced production level of specific volatiles.
Pseudoxylaria antagonizes Termitomyces growth and metabolizes fungal biomassThe growth behavior of Pseudoxylaria isolates was also analyzed in co-culture assays with Termitomyces. As expected from prior studies, both fungi showed reduced growth when co-cultured on agar plates, often causing the formation of zones of inhibition (ZOI) between the fungal colonies (Fig. 6A–D, Table S11–S14) [7]. When fungus-fungus co-cultures were maintained for longer than two weeks on agar plates, Pseudoxylaria started to overcome the ZOI and overgrew Termitomyces via the extension of aerial mycelium. The observation was even more pronounced when co-cultures were performed on wood-rice medium (WRM), where Pseudoxylaria remained the only visible fungus after two weeks.Fig. 6: Co-cultivation of Pseudoxylaria sp. X170LB and Termitomyces sp. T112 and results of isotope fractionation experiments.Representative pictures of fungal growth and co-cultivation of A Termitomyces sp. T112, B Pseudoxylaria sp. X170LB, C co-culture of Pseudoxylaria sp. X802 and Termitomyces sp. T153 exhibiting a ZOI, in which X802 overgrowths T153 in proximity to the interaction zone (red arrow), and D Pseudoxylaria sp. X802 growing on the surface of a living Termitomyces sp. T153 culture. E, F Shown is the relative change in the carbon isotope pattern (δ13C values, ± standard deviation, with n = 3) of lipid and carbohydrate fractions isolated from fungal biomass of Termitomyces sp. T112, Pseudoxylaria sp. X170LB, and Pseudoxylaria sp. X170LB cultivated on vegetative Termitomyces sp. T112 biomass (T112ǂ), or on lyophilized Termitomyces sp. T112 biomass (T112). Fungal strains were grown on E medium with natural 13C abundance and F medium artificially enriched in 13C content.Full size imageTo verify whether Pseudoxylaria consumes Termitomyces or even partially degrades specific metabolites present within the fungal biomass, we pursued stable isotope fingerprinting commonly used to analyse trophic relations [38, 39]. This diagnostic method relies on measurable changes in the bulk stable isotope composition, because biosynthetic enzymes preferentially convert lighter metabolites enriched in 12C compared to their heavier 13C-enriched congeners. This intrinsic kinetic isotope effect results in an overall change in the 13C/12C ratio of the respective educts and products, in particular in biomarkers such as phospholipid fatty acids, carbohydrates, and amino acids. Using this isotope enrichment effect, we determined the natural trophic isotope fractionation of 13C in lipids and carbohydrates produced by Termitomyces sp. T112 and Pseudoxylaria sp. X170LB. For clearer differentiation, both fungi were cultivated on PDA medium containing naturally abundant 13C/12C, Fig. 6E) and on PDA medium enriched with 13C-glucose (Fig. 6F). Lipids and carbohydrates were isolated from mycelium harvested after 21 days (Fig. 6E, Table S15).Analysis of fungal carbohydrate and lipid-rich metabolite fractions by Elemental Analysis-Isotope Ratio Mass Spectrometry (EA-IRMS) [40, 41] uncovered that under normal growth conditions (full medium), Termitomyces sp. T112 and Pseudoxylaria sp. X170LB showed only a slight negative trophic fractionation of stable carbon isotopes (δ13C/12C ratio (expressed as δ13C values [‰]), Fig. 6F) within the carbohydrate fractions (T112: −1.2 ‰; for X170LB: −1.3 ‰), and expectedly a stronger depletion in the lipid fraction (T112: −6.7 ‰, and less pronounced for X170LB: −3.1 ‰). To determine if Pseudoxylaria metabolizes Termitomyces biomass, the isotope pattern of metabolites derived from Pseudoxylaria thriving on living biomass of Termitomyces (T112ǂ) was analysed next. Here, an overall positive carbon isotope (13C/12C) fractionation by approximately +0.6 ‰ relative to the control medium was detectable, while the δ13C values of lipids remained largely unchanged (Fig. 6F, Table S15). These results suggested that Pseudoxylaria might pursue a preferential uptake of Termitomyces-derived carbohydrates.In a last experiment, Pseudoxylaria was grown on lyophilized (dead) Termitomyces biomass (T112) as sole food source. In this experiment, the isotope fingerprint showed converging δ13C values of −1.9 ‰ (relative to the media) for both carbohydrate and lipid fractions, which indicated that Pseudoxylaria is able to simultaneously metabolize and cycle carbohydrates as well as lipids resulting in the equilibration of isotopic levels between carbohydrates and lipids. Thus, it was concluded that in nature, Pseudoxylaria likely harvests nutrients firstly from vegetative Termitomyces, and then—if possible—subsequently degrades dying or dead mycelium.
Pseudoxylaria produces antimicrobial secondary metabolitesBased on the observation that Pseudoxylaria antagonizes growth of Termitomyces, we questioned if the formation of a ZOI might be caused by the secretion of Pseudoxylaria-derived antimicrobial metabolites [26, 42]. Thus, we performed an ESI(+)-HRMS/MS based metabolic survey using the web-based platform “Global Natural Product Social Molecular Networking” (GNPS) [43] to correlate the encoded biosynthetic repertoire of Pseudoxylaria with secreted metabolites.A partial similar metabolic repertoire across the six analyzed strains was detectable and allowed us to match some of the detectable chemical features and previously isolated metabolites to the predicted shared BGCs, such as antifungal and histone deacetylase inhibitory xylacremolides (Xyl; X187/Mn132) [32, 33], pseudoxylaramides (Psa; X187/Mn132) [32], antibacterial pseudoxylallemycins (Psm; X802/OD126) [18], xylasporin/cytosporins (Xsp; X802/OD126/X187/Mn132) [36], and cytotoxic cytochalasins (X802/OD126) (Fig. 7A and B) [18].Fig. 7: Comparative metabolomic analysis of six Pseudoxylaria strains (OD126 (red), Mn132 (orange), X170 (black), X187 (green), X3.2 (yellow) and X802 (blue)).A Overview of the GNPS network. Identified metabolite clusters xylacremolides (Xyl; X187/Mn132) [32, 33], pseudoxylaramides [32] (Psa; X187/Mn132), pseudoxylallemycins (Psm; X802/OD126) [18], xylasporin/cytosporins (Xsp; X802/OD126/X187/Mn132) and cytochalasins (X802/OD126) [18]. B xylasporin/cytosporin-related cluster formed by nodes from X802 (blue), OD126 (red), X187 (green) and Mn132 (orange). C Chemical structures of natural products isolated from Pseudoxylaria species and related compounds. Red box highlights proposed structures of isolated xylasporin G and I in this study.Full size imageA cluster that contained MS2 signals of molecular ions assigned to the cytosporin/xylasporin family, which was shared by at least four strains, caught our attention as a certain degree of structural diversity of xylasporin/cytosporin family was predicted from the comparison of their respective BGCs. The assigned nodes of this GNPS cluster split into two subclusters with only very little overlap between both regions. Analysis of the mass fragment shifts suggested that both subclusters belong to two different families of xylasporin/cytosporin congeners (Figure S9). To verify these deductions, we pursued an MS-guided purification of xylasporin/cytosporins from chemical extracts of Pseudoxylaria sp. X187, which yielded xylasporin G (3.23 mg, pale-yellow solid) and xylasporin I (1.75 mg, pale-yellow solid). The sum formulas of xylasporin G and xylasporin I were determined to be C17H22O5 (calcd. for [M + H]+ C17H23O5+ = 307.1540, found 307.15347, −1.726 ppm) and C17H24O5 (calcd. for [M + H]+ C17H25O5+ = 309.1697, found 309.1691, −1.68 ppm) by ESI-(+)-HRMS and were predicted to have six degrees of unsaturation (Fig. 7B, Figure S10, Table S16-S17). Planar structures were deduced by comparative 1D and 2D NMR analyses, which revealed the presence of an unsaturated polyketide chain that matched the unsaturation degree and the anticipated structural variation from cytosporins (Fig. 7C, Figure S11-S25).To evaluate if Pseudoxylaria-derived culture extracts and produced natural products (e.g., cytochalasins) are responsible for the observed antimicrobial activity, standardized antimicrobial activity assays were performed (Table S17, S18 and Figure S26). As neither culture extracts nor single compounds exhibited significant antimicrobial activity, they could not be held fully accountable for the antagonistic behavior in co-cultures. Thus, we hypothesized that the observed ZOI might be caused by yet unknown effects like nutrient depletion or bioactive enzymes.
Pseudoxylaria has a negative impact on the fitness of insect larvaeDue to the production of structurally diverse and weakly antimicrobial secondary metabolites, we questioned if mycelium of Pseudoxylaria exhibits intrinsic insecticidal or other insect-detrimental activities, which could discourage or ward off grooming behavior of termite workers. Due to the technical challenges associated with behavioral studies of termites, we evaluated instead the effect of Pseudoxylaria biomass on Spodoptera littoralis, a well-established insect model system and a destructive agricultural lepidopterous pest [44, 45]. When S. littoralis larvae were fed with mycelium-covered agar plugs of Pseudoxylaria sp. X802, a clear decrease of the relative growth rate (RGR) and decline in survival was observed (Fig. 8: treatment D (green), Table S19, S20) compared to feeding with untreated agar plugs (treatment A (black)). In comparison, when larvae were fed with agar plugs covered with the fungal mutualist Termitomyces sp. T153 (treatment B (blue)) an increased growth rate of larvae was observed.Fig. 8: Effect of Termitomyces sp. T153 and Pseudoxylaria sp. X802 mycelia on the relative growth rate and survival of S. littoralis larvae.Insects were fed with either A PDA, B PDA agar plug covered with vegetative Termitomyces sp. T153, C PDA agar plug from which vegetative Termitomyces sp. T153 was removed prior to feeding, D PDA agar plug covered with vegetative Pseudoxylaria sp. X802 mycelium, and E PDA agar plug from which vegetative Pseudoxylaria sp. X802 mycelium was removed prior to feeding. All experiments were performed with 25 replicates per treatment, a duration of 10 days, and larval weights and survival rates were recorded every day. Statistical significances were determined using ANOVA on ranks (p  More