# Plant pathogen life-history traits and adaptation to environmental constraints

Christophe Le May, Josselin Montarry, Cindy E. Morris, Omer Frenkel, Virginie Ravigné

► To cite this version:

Christophe Le May, Omer Frenkel, Virginie Ravigné (Dir.). Plant pathogen life-history traits and adaptation to environmental constraints. Frontière Media, 126 p., 2020, Frontiers Research Topics, 978-2-88963-530-6. 10.3389/978-2-88963-530-6. hal-05177356

HAL Id: hal-05177356

https://hal.science/hal-05177356v1

Submitted on 22 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
PLANT PATHOGEN LIFE-HISTORY
TRAITS AND ADAPTATION TO
ENVIRONMENTAL CONSTRAINTS

EDITED BY: Christophe Le May, Josselin Montarry, Cindy E. Morris,
Omer Frenkel and Virginie Ravigné

PUBLISHED IN: Frontiers in Plant Science and Frontiers in Microbiology

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

frontiers

Frontiers eBook Copyright Statement

The copyright in the text of
individual articles in this eBook is the
property of their respective authors
or their respective institutions or
funders. The copyright in graphics
and images within each article may
be subject to copyright of other
parties. In both cases this is subject
to a license granted to Frontiers.

The compilation of articles
constituting this eBook is the
property of Frontiers.

Each article within this eBook, and
the eBook itself, are published under
the most recent version of the
Creative Commons CC-BY licence.
The version current at the date of
publication of this eBook is
CC-BY 4.0. If the CC-BY licence is
updated, the licence granted by
Frontiers is automatically updated to
the new version.

When exercising any right under the
CC-BY licence, Frontiers must be
attributed as the original publisher
of the article or eBook, as
applicable.

Authors have the responsibility of
ensuring that any graphics or other
materials which are the property of
others may be included in the
CC-BY licence, but this should be
checked before relying on the
CC-BY licence to reproduce those
materials. Any copyright notices
relating to those materials must be
complied with.

Copyright and source
acknowledgement notices may not
be removed and must be displayed
in any copy, derivative work or
partial copy which includes the
elements in question.

All copyright, and all rights therein,
are protected by national and
international copyright laws. The
above represents a summary only.

For further information please read
Frontiers' Conditions for Website
Use and Copyright Statement, and
the applicable CC-BY licence.

ISSN 1664-8714

ISBN 978-2-88963-530-6

DOI 10.3389/978-2-88963-530-6

About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view.

By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org
PLANT PATHOGEN LIFE-HISTORY
TRAITS AND ADAPTATION TO
ENVIRONMENTAL CONSTRAINTS

Topic Editors:

**Christophe Le May**, Agrocampus Ouest, France

**Josselin Montarry**, INRAe Centre Bretagne-Normandie, France

**Cindy E. Morris**, INRAe Centre Provence-Alpes-Côte d'Azur, France

**Omer Frenkel**, Agricultural Research Organization (ARO), Israel

**Virginie Ravigné**, Centre de Coopération Internationale en Recherche
Agronomique pour le Développement (CIRAD), France

**Citation:** Le May, C., Montarry, J., Morris, C. E., Frenkel, O., Ravigné, V., eds. (2020). Plant Pathogen Life-History Traits and Adaptation to Environmental Constraints. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-530-6
# Table of Contents

**04 Editorial: Plant Pathogen Life-History Traits and Adaptation to Environmental Constraints**
Christophe Le May, Josselin Montarry, Cindy E. Morris, Omer Frenkel and Virginie Ravigné

**08 Genome Scans Reveal Homogenization and Local Adaptations in Populations of the Soybean Cyst Nematode**
Anne-Frédérique Gendron St-Marseille, Etienne Lord, Pierre-Yves Véronneau, Jacques Brodeur and Benjamin Mimee

**22 Real-Time PCR for Diagnosing and Quantifying Co-infection by Two Globally Distributed Fungal Pathogens of Wheat**
Araz S. Abdullah, Chala Turo, Caroline S. Moffat, Francisco J. Lopez-Ruiz, Mark R. Gibberd, John Hamblin and Ayalsew Zerihun

**32 Defense Compounds Rather Than Nutrient Availability Shape Aggressiveness Trait Variation Along a Leaf Maturity Gradient in a Biotrophic Plant Pathogen**
Agathe Maupetit, Romain Larbat, Michaël Pernaci, Axelle Andrieux, Cécile Guinet, Anne-Laure Boutigny, Bénédicte Fabre, Pascal Frey and Fabien Halkett

**45 Adaptation of a Fungal Pathogen to Host Quantitative Resistance**
Lise Frézal, Guy Jacqua and Claire Neema

**57 Populations of the Beet Cyst Nematode Heterodera schachtii Exhibit Strong Differences in Their Life-History Traits Across Changing Thermal Conditions**
Sylvain Fournet, Lucile Pellan, Catherine Porte, Christophe Piriou, Eric Grenier and Josselin Montarry

**66 Role of Temperature and Coinfection in Mediating Pathogen Life-History Traits**
Elise Vaumourin and Anna-Liisa Laine

**74 Striking Similarities Between Botrytis cinerea From Non-agricultural and From Agricultural Habitats**
Marc Bardin, Christel Leyronas, Claire Troulet and Cindy E. Morris

**85 Pathogenicity Traits Correlate With the Susceptible Vitis vinifera Leaf Physiology Transition in the Biotroph Fungus Erysiphe necator: An Adaptation to Plant Ontogenic Resistance**
Agnes Calonnec, Jerome Jolivet, Philippe Vivin and Sylvain Schnee

**102 A Trade-Off Between Sporangia Size and Number Exists in the Potato Late Blight Pathogen Phytophthora infestans, and is not Altered by Biotic and Abiotic Factors**
Nicolas Mariette, Alexander Kröner, Romain Mabon, Josselin Montarry, Bruno Marquer, Roselyne Corbière, Annabelle Androdias and Didier Andrivon

**109 Responses to Hydric Stress in the Seed-Borne Necrotrophic Fungus Alternaria brassicicola**
Guillaume Quang N'Guyen, Roxane Raulo, Muriel Marchi, Carlos Agustí-Brisach, Beatrice Iacomi, Sandra Pelletier, Jean-Pierre Renou, Nelly Bataillé-Simoneau, Claire Campion, Franck Bastide, Bruno Hamon, Chloé Mouchès, Benoit Porcheron, Remi Lemoine, Anthony Kwasiborski, Philippe Simoneau and Thomas Guillemette
<figure><img src="image_4.png" /></figure>

Check for
updates

# Editorial: Plant Pathogen Life-History Traits and Adaptation to Environmental Constraints

Christophe Le May¹*, Josselin Montarry¹, Cindy E. Morris², Omer Frenkel³ and Virginie Ravigné⁴

¹IGEPP, INRA, Agrocampus-Ouest, Université de Rennes 1, Le Rheu, France, ²INRA, UR Pathologie Végétale, Avignon, France, ³Agricultural Research Organization, Volcani Center, Rishon LeZion, Israel, ⁴BIOS Departement, PVBMT Laboratory, CIRAD, Montpellier, France

**Keywords:** coinfection, epidemic, host resistance, hydric stress, leaf maturity, non-agricultural habitats, survival, temperature

## Editorial on the Research Topic

### Plant Pathogen Life-History Traits and Adaptation to Environmental Constraints

Life-history is a key concept in evolutionary biology and ecology. It corresponds to the narrative of the various events punctuating the existence of an organism from his birth to his death (Begon et al., 2006; Michalakis et al., 2016). Throughout life, living organisms acquire resources that they actively find or extract from their environment and then allocate to different functions: development, survival, dispersal and reproduction (Roff, 1992; Stearns, 1992; Michalakis, 2009). Traits involved in the timing and amplitude of these allocation dilemmas are defined as life history traits. Life-history traits are often determinant for individuals (Kingsolver and Pfennig, 2007). They influence spatial and temporal disease dynamics, and thus the genetic diversity and structure of pathogen populations (Barrett et al., 2008; Michalakis et al., 2016). These determinants are involved in the ability of pathogens to adapt to varying ecological factors including changes in the biotic and abiotic effect, but also direct or indirect interactions with other strains or species of pathogens co-occurring on the same host (Michalakis et al., 2016; Tollenaere et al., 2016). For plant pathogens, these traits can be grouped into two categories: those involved in the epidemic phase, onto or into the host, and those related to the survival phase, often outside the host. Trade-off can occurred between these two phases, leading to consequences on epidemic dynamics, evolution, and speciation (Pariaud et al., 2009; Hamelin et al., 2011). Understanding processes maintaining variation in plant pathogen life-history traits is a central question in evolutionary ecology and a major challenge for the design of disease control strategies (Galvani, 2003; Grenfell et al., 2004).

The aim of this Research Topic was to collect research papers with vast perspective on plant pathogen life-history traits and adaptation to environmental constraints. Life-history traits of both epidemic and survival phases were considered, and tested environmental constraints included abiotic (temperature, precipitation, hydric stress) and biotic (habitat of origin, cultivar, maturity of leaves, coinfection) factors. Among the 10 accepted papers, three types of plant pathogens were considered: fungi (seven papers), nematodes (two papers), and oomycetes (one paper).

Regarding the epidemic phase, four main topics were covered by the research papers accepted in this special issue: i) host adaptation, ii) impact of host's physiological status, iii) adaptation to climatic conditions, and iv) effects of coinfection.

## OPEN ACCESS

**Edited and reviewed by:**
Brigitte Mauch-Mani,
Université de Neuchâtel,
Switzerland

*Correspondence:
Christophe Le May
lemay@agrocampus-ouest.fr

**Specialty section:**
This article was submitted to
Plant Microbe Interactions,
a section of the journal
Frontiers in Plant Science

Received: 07 November 2019

Accepted: 09 December 2019

Published: 24 January 2020

**Citation:**
Le May C, Montarry J, Morris CE, Frenkel O and Ravigné V (2020)
*Editorial: Plant Pathogen Life-History Traits and Adaptation to Environmental Constraints.*
*Front. Plant Sci.* 10:1730.
doi: 10.3389/fpls.2019.01730
How host adaptation affects life history traits evolution was particularly analyzed by studying the effect of resistant cultivars onto several life-history traits. In the ascomycete *Colletotrichum gloeosporioides* the impact of host quantitative resistance on pathogen evolution was studied by evaluating the diversity of populations at the field scale at both neutral markers and pathogenic traits (Frezal et al.). Population genetic structure revealed a significant influence of clonal reproduction in *C. gloeosporioides* evolution and a low migration rate between fields. Results of cross-inoculation tests showed that aggressiveness of the fungal clones seemed to have evolved through an accumulation of components specific to each water yam cultivar, suggesting an adaptation to their host cultivar. Despite the remaining marks of adaptation to the former widely cultivated host, adaptation to current cultivars was clearly depicted. Such pattern of local adaptation to the dominant cultivar was previously highlighted in the potato late blight pathogen *Phytophthora infestans* (Andrivon et al., 2007). In their study, Mariette et al. explored whether the now classical negative relationship between offspring size and number occurred in this oomycete, *P. infestans*, and whether the trade-off was impacted by potato cultivar or host of origin of the pathogen (tomato and potato). They confirmed the existence of a trade-off and showed that it was not affected by any of these biotic factors. The observed polyphenism for these traits in *P. infestans* populations will favor the coexistence of distinct reproductive strategies.

Over the last years, intra-host dynamics has been identified as a major component of epidemiological and evolutionary dynamics of pathogens (Alizon et al., 2009). In plants it is known that host tissue characteristics vary a lot with tissue age. In particular, according the age of tissues, nutrient availability as well as defense levels may vary significantly, with opposing effects on the success of infection (Kus et al., 2002; Al-Naimi et al., 2005). Ontogenic resistance against powdery mildew (*Erysiphe necator*) on *Vitis vinifera* was studied by Calonnec et al. They showed that the three pathogenic traits studied (infection efficiency, sporulation, and mycelium growth) were affected, and their variation was strongly correlated with leaf age. Sporulation was more closely correlated with variations in sugar contents and the infection efficiency with leaf water content, suggesting that ontogenic resistance on grapevine leaves seems to be an immutable physiological process that *E. necator* is able to circumvent by restricting its development to sink tissue. Similarly, Maupetit et al. tested how nutrient availability and defense level impacted the development of *Melampsora larici-populina* (infection efficiency, latent period, uredinia size, mycelium quantity, sporulation rate, sporulation capacity, and spore volume). They showed that *M. larici-populina* was more aggressive on more mature leaves as indicated by wider uredinia and a higher sporulation rate. In contrast, phenolic contents (flavonols, hydroxycinnamic acid esters, and salicinoids) were negatively correlated with uredinia size and sporulation rate suggesting that pathogen's fitness appeared to be more constrained by the constitutive plant defense level than limited by nutrient availability, as evident in the decrease in sporulation.

Like other organisms, pathogens need a particular range of temperature and humidity in order to invade and develop on their host plant (Agrios, 2005). These two parameters particularly can modify spore germination, host colonization, spore production and dispersal of fungal pathogens (Fitt et al., 2006). Decades of research have generated considerable knowledge and greater understanding of the seasonal effects of temperature, rainfall, and humidity on diseases affecting major food crops (Garrett et al., 2011; Pautasso et al., 2011; Donatelli et al., 2017). In this special issue, two main parameters were studied, change in temperature and effect of water stress. Regarding the effect of temperature, Vaumourin and Laine developed a study on *Podosphaera plantaginis* and showed that temperature had a significant effect on all measured life-history traits, suggesting that the effect of temperature on life-history traits was both direct as well as mediated through a genotype-by-temperature interaction. Similarly, Mariette et al. showed on *P. infestans* that the negative relationship between offspring size and number was maintained for all the tested temperatures. Using cyst nematodes, other studies showed an effect of temperature on different life-history traits. Indeed, Fournet et al. showed a strong temperature effect on life-history traits of the beet cyst nematode *Heterodera schachtti*. While nematode multiplication was not differentially affected by temperatures, as favorable conditions for the host are also favorable for the parasite, the effect of temperature on hatching depended on the origin of populations, separating southern from northern European populations. Using a genome scan approach, Gendron St-Marseille et al. evaluated the adaptive potential of the soybean cyst nematode *Heterodera glycines*. The genotyping by sequencing of 64 *H. glycines* populations, allowed identifying 15 loci under selection for climatic or geographic co-variables. Lastly, regarding the effect of water stress, a transcriptomic and metabolomic approach was developed on the *Arabidopsis thaliana/Alternaria brassicicola* pathosystem, to study the level of susceptibility of the fungus to water stress and its impact on its seed transmission ability (N'Guyen et al.). These approaches led to the identification of specific proteins (hydrolipin) implicated in the tolerance toward water stress, and seed transmission, but not in the aggressiveness of the pathogen.

The ability of a pathogen to establish and grow on its host may be drastically altered by simultaneous infection by other pathogens (other strains or species of pathogens co-occurring on the same host). Life-history allocation may change under coinfection, affecting the evolutionary potential and epidemiological dynamics of pathogens (Tollenaere et al., 2016; Laine and Mäkinen, 2018). In their study, Vaumourin and Laine found that coinfection only modified the number of sexual resting structures produced, but did not changed the production of asexual forms. Coinfection studies, assessment of fungi sharing common characteristics, and host species, create a challenge for conventional disease diagnosis and subsequent management strategies. Indeed, lack of powerful tools to simultaneously distinguish and quantify these types of pathogens constitutes a major limitation. In their study, Abdullah et al. developed a duplex real-time PCR assay for
quantifying of co-infecting wheat pathogens, *Pyrenophora tritici-repentis* and *Parastagonospora nodorum*. The utility of the method was demonstrated using field samples of a cultivar sensitive to both pathogens. While visual and culture diagnosis suggested the presence of only one of the pathogen species, the assay revealed not only presence of both co-infecting pathogens (hence enabling asymptomatic detection) but also allowed quantification of relative abundances of the pathogens as a function of disease severity.

Survival stage outside the host has often been much overlooked due to the difficulty of its study. It is a key stage in pathogen's life cycle. Pathogens are all dependent on the life cycle, where despite the existence of dormant forms, evolutionary processes may occur. Thus, according to the status of their host (perennial or annual host plant), pathogens have to develop a wide array of strategies to overcome the alternation of presence/absence (and/or modification of physiological activity) of their host (Agrios, 2005). This survival phase is consequently important as it may affect the population size of the pathogen and lead to a bottleneck process. In this, special issue, research developed by Bardin et al. have investigated the ubiquity of the broad host range necrotrophic fungus *Botrytis cinerea*, outside of agricultural settings and have determined whether the populations in these natural habitats can be distinguished phenotypically and phylogenetically from populations isolated from diseased crops. Their results showed that *B. cinerea* strains sampled on different non-agricultural substrates were genetically and phenotypically similar to

strains sampled in agricultural substrates. These results suggest that highly diverse populations of this plant pathogen persist outside of agriculture in association with substrates other than cultivated plants and that this component of their life-history is compatible with its capability to maintain its potential as plant pathogen. The survival phase was also considered in Vaumourin and Laine, as coinfection affects the number of sexual resting structures produced, and in Fournet et al., as they considered also the hatching after storage at different temperatures, simulating survival conditions during the inter-cropping period.

Considering plant pathogen traits in the context of life history trait evolution is a promising research avenue, still in its infancy. We hope that this research topic will help foster this line of research and provide a valuable resource for researchers working in the field of evolutionary ecology of plant pathogen populations.

## AUTHOR CONTRIBUTIONS

CL, JM, CM, OF, and VR have made substantial direct contributions to the work, and approved it for publication.

## ACKNOWLEDGMENTS

We thank all the authors involved in this research topic for their excellent and varied contributions.

## REFERENCES

Agrios, G. N. (2005). *Plant Pathology- 5th ed* (United States: Elsevier Academic Press). 2005.

Alizon, S., Hurford, A., Mideo, N., and van Baalen, M. (2009). Virulence evolution and the trade-off hypothesis: history, current state of affairs and the future. *J. Evol. Biol.* 22, 245–259. doi: 10.1111/j.1420-9101.2008.01658.x

Al-Naimi, F. A., Garrett, K. A., and Bockus, W. W. (2005). Competition, facilitation, and niche differentiation in two foliar pathogens. *Oecologia* 143, 449–457. doi: 10.1007/s00442-004-1814-x

Andrivon, D., Pilet, F., Montarry, J., Hafidi, M., Corbière, R., Achbani, E. H., et al. (2007). Adaptation of *Phytophthora infestans* to partial resistance in potato: evidence from French and Moroccan populations. *Phytopathol* 97, 338–343. doi: 10.1094/PHYTO-97-3-0338

Barrett, L. G., Thrall, P. H., Burdon, J. J., and Linde, C. C. (2008). Life history determines genetic structure and evolutionary potential of host–parasite interactions. *Trends Ecol. Evol.* 23, 678–685. doi: 10.1016/j.tree.2008.06.017

Begon, M., Townsend, C. R., and Harper, J. L. (2006). *Ecology: from individuals to ecosystems* (Malden, MA: Blackwell Pub.).

Donatelli, M., Magarey, R. D., Bregaglio, S., Willocquet, L., Whish, J. P. M., and Savary, S. (2017). Modelling the impacts of pests and diseases on agricultural systems. *Agric. Syst.* 155, 213–224. doi: 10.1016/j.agsy. 2017.01.019

Fitt, B. D. L., Huang, Y.-J., van den, Bosch F, and JS, West (2006). Coexistence of related pathogen species on arable crops in space and time. *Ann. Rev. Phytopathol.* 44, 163–182. doi: 10.1146/annurev.phyto.44.070505.143417

Galvani, A. P. (2003). Epidemiology meets evolutionary ecology. *Trends Ecol. Evol.* 18 (3), 132–139.

Garrett, K. A., Nita, M., De Wolf, E. D., Esker, P. D., Gomez-Montano, L., and Sparks, A. H. (2011). “Plant pathogens as indicators of climate change,” in

*Climate Change (2nd edn): Observed Impacts on Planet Earth*. Ed. Trevor M. Letcher (Amsterdam, The Netherlands: Elsevier BV), 325–338.

Grenfell, B. T., Pybus, O. G., Gog, J. R., Wood, J. L. N., Daly, J. N., Mumford, J. A., and Holmes, E. C. (2004). Unifying the epidemiological and evolutionary dynamics of pathogens. *Science* 303 (5656), 327–332.

Hamelin, F., Castel, M., Poggi, S., Andrivon, D., and Mailleret, L. (2011). Seasonality and the evolutionary divergence of plant parasites. *Ecology* 92 (12), 2159–2166. doi: 10.1890/10-2442.1

Kingsolver, J. G., and Pfennig, D. W. (2007). Patterns and power of phenotypic selection in nature. *BioScience* 57 (7), 561–572. doi: 10.1641/B570706

Kus, J. V., Zaton, K., Sarkar, R., and Cameron, R. K. (2002). Age-related resistance in Arabidopsis is a developmentally regulated Ddfense response to *Pseudomonas syringae*. *Plant Cell* 14, 479–490. doi: 10.1105/tpc.010481

Laine, A.-L., and Mäkinen, H. (2018). Life-history correlations change under coinfection leading to higher pathogen load. *Evol. Lett.* 2, 126–133. doi: 10.1002/evl3.48

Michalakis, Y., Charmantier, A., Gaillard, J.-M., Sorci, G., Tully, T., and Ronce, O. (2016). “Evolution des traits d’histoire de vie,” in *Biologie évolutive* (Louvain-la-Neuve: De Boeck Supérieur), 373–422.

Michalakis, Y. (2009). “Parasitism and the evolution of life-history traits,” in *Ecology and evolution of parasitism* (Oxford: Oxford University Press), 19–30.

Pariaud, B., Ravigne, V., Halkett, F., Goyeau, H., Carlier, J., and Lannou, C. (2009). Aggressiveness and its role in the adaptation of plant pathogens. *Plant Pathol.* 58, 409–424. doi: 10.1111/j.1365-3059.2009.02039.x

Pautasso, M., Doring, T. F., Garbelotto, M., Pellis, L., and Jeger, M. J. (2011). Impacts of climate change on plant diseases—opinions and trends. *Eur. J. Plant Pathol.* 133, 295–313. doi: 10.1007/s10658-012-9936-1

Roff, D. A. (1992). *The evolution of life histories: theory and analysis* (New York: Chapman and Hall).

Stearns, S. C. (1992). *The evolution of life histories* (Oxford: Oxford University Press).
Tollenaere, C., Susi, H., and Laine, A.-L. (2016). Evolutionary and epidemiological implications of multiple infection in plants. *Trends Plant Sci.* 21, 80–90. doi: 10.1016/j.tplants.2015.10.014

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Le May, Montarry, Morris, Frenkel and Ravigné. This is an open-
access article distributed under the terms of the Creative Commons Attribution
License (CC BY). The use, distribution or reproduction in other forums is permitted,
provided the original author(s) and the copyright owner(s) are credited and that the
original publication in this journal is cited, in accordance with accepted academic
practice. No use, distribution or reproduction is permitted which does not comply with
these terms.
<figure><img src="image_5.png" /></figure>

Check for
updates

# Genome Scans Reveal Homogenization and Local Adaptations in Populations of the Soybean Cyst Nematode

Anne-Frédérique Gendron St-Marseille¹,², Etienne Lord¹, Pierre-Yves Véronneau¹, Jacques Brodeur² and Benjamin Mimee¹*

¹ Saint-Jean-sur-Richelieu Research and Development Centre, Agriculture and Agri-Food Canada, Saint-Jean-sur-Richelieu, QC, Canada, ² Institut de Recherche en Biologie Végétale (IRBV), Université de Montréal, Montréal, QC, Canada

Determining the adaptive potential of alien invasive species in a new environment is a key concern for risk assessment. As climate change is affecting local climatic conditions, widespread modifications in species distribution are expected. Therefore, the genetic mechanisms underlying local adaptations must be understood in order to predict future species distribution. The soybean cyst nematode (SCN), *Heterodera glycines* Ichinohe, is a major pathogen of soybean that was accidentally introduced in most soybean-producing countries. In this study, we explored patterns of genetic exchange between North American populations of SCN and the effect of isolation by geographical distance. Genotyping-by-sequencing was used to sequence and compare 64 SCN populations from the United States and Canada. At large scale, only a weak correlation was found between genetic distance (Wright's fixation index, $F_{ST}$) and geographic distance, but local effects were strong in recently infested states. Our results also showed a high level of genetic differentiation within some populations, allowing them to adapt to new environments and become established in new soybean-producing areas. Bayesian genome scan methods identified 15 loci under selection for climatic or geographic co-variables. Among these loci, two non-synonymous mutations were detected in SMAD-4 (mothers against decapentaplegic homolog 4) and DOP-3 (dopamine receptor 3). High-impact variants linked to these loci by genetic hitchhiking were also highlighted as putatively involved in local adaptation of SCN populations to new environments. Overall, it appears that strong selective pressure by resistant cultivars is causing a large scale homogenization with virulent populations.

**Keywords:** Bayesian outlier detection, genetic diversity, genotyping-by-sequencing, *Heterodera glycines*, isolation by distance

## INTRODUCTION

The introduction of an organism into a new environment can have unpredictable detrimental consequences, including public health problems, losses in biodiversity and ecosystem services or crop yield losses due to exotic weeds, insects and pathogens, altogether resulting in significant economic impact (Pimentel et al., 2005). Unfortunately, the steady increase in international trade

**OPEN ACCESS**

**Edited by:**

Christophe Le May,
Agrocampus Ouest, France

**Reviewed by:**

Yessica Rico,
Instituto de Ecología (INECOL), Mexico
Sebastian Eves-van Den Akker,
University of Cambridge,
United Kingdom

**\*Correspondence:**

Benjamin Mimee
benjamin.mimee@agr.gc.ca

**Specialty section:**

This article was submitted to
Plant Microbe Interactions,
a section of the journal
Frontiers in Plant Science

**Received:** 21 March 2018

**Accepted:** 18 June 2018

**Published:** 17 July 2018

**Citation:**

Gendron St-Marseille A-F, Lord E, Véronneau P-Y, Brodeur J and Mimee B (2018) Genome Scans Reveal Homogenization and Local Adaptations in Populations of the Soybean Cyst Nematode. *Front. Plant Sci.* 9:987. doi: 10.3389/fpls.2018.00987
facilitates the movement and introduction of new invasive species (Hulme, 2009). In addition, climate change is altering environmental conditions and could change the species' distribution range or favor their establishment in previously unsuitable habitats (Early and Sax, 2014). It is therefore imperative to carry out specific risk assessment in order to target species to be controlled. Alien invasive species, by definition, did not evolve in the biogeographic habitat in which they are introduced. Consequently, they are often poorly adapted to their new environment. Generally, the most successful invaders will have a high potential for rapid adaptation through phenotypic plasticity or microevolution (Novak, 2007). Understanding the genetic mechanisms of local adaptation is therefore critical to predict future species distribution.

Plant-parasitic nematodes are microscopic worms that reduce global annual food production by 12.3% and cause more than US$157 billion in economic losses worldwide (Hassan et al., 2013). In North America, one of the most damaging species is the soybean cyst nematode (SCN), *Heterodera glycines* Ichinohe. Since its first detection in 1954 in Hanover County, North Carolina, SCN has been reported in almost every soybean-producing county in the United States (Davis and Tylka, 2000), as well as in southwestern Ontario, Canada (Anderson et al., 1988). During the 2000s, SCN colonized several new northern localities in the United States (North Dakota in 2003; Mathew et al., 2015) and Eastern Canada (northeastern Ontario in 2007 and Québec in 2013; Mimee et al., 2014). A few studies have investigated SCN dispersal, which seems to follow a northward and eastward pathway in North America (Tylka and Marett, 2017). This dispersion pattern correlates with the expansion of soybean cultivation following the introduction of new soybean varieties with shorter maturity periods and improved tolerance to drought and cold (Shurtleff and Aoyagi, 2010; Yu, 2011).

It is expected that SCN would survive and multiply throughout the current North American soybean-growing area and complete at least two generations at its northern limit (Gendron St-Marseille, 2013). Climate warming should also favor the establishment of SCN at higher latitudes and increase the number of generations per year in most regions. However, these predictions are based strictly on temperature requirements and do not account for the intrinsic capacity of SCN to adapt to new environmental conditions. Genetic variations within a population reflect its evolutionary potential and result from four evolutionary forces that affect individual fitness: mutation, gene flow, selection, and genetic drift (Eizaguirre and Baltazar-Soares, 2014). For most organisms, including SCN, the relative weights of these forces can differ significantly. Mutations are rare events that should not contribute significantly to SCN adaptations in a short time frame. Gene flow depends on the dispersal ability of an organism, which for SCN is achieved mainly by means of human activities at the short spatial scale (Kristjansson, 2010). In addition, wind and flooding can carry SCN cysts over very long distances, and contribute to its dispersal at the regional and continental scales. Many different selection pressures can shape the genetic structure of SCN populations, but host plant is probably the strongest selection factor. The ability of SCN to reproduce on a given soybean genotype differs greatly depending

on its resistance genes and the nematode's virulence profile (HG type) (Colgrove and Niblack, 2008; Niblack et al., 2008). Thus, management decisions by growers (for example, the systematic use of resistant cultivars) can result in a strong selection pressure. Finally, the influence of genetic drift will also depend on pest control strategies, because they contribute to dictate population size, although it was shown for cyst nematodes that genetic diversity can be very high within a single cyst (Green et al., 1970). Each nematode female can mate with several males and lay hundreds of eggs that can survive for at least a decade in the soil (Slack et al., 1972). Thus, even if the diversity appears reduced due to genetic drift under strong selection by resistant cultivars, most alleles probably persist for several years in the population in a "dormant" state.

In other cyst nematode species, genetic diversity at the population level has been studied by means of several techniques, including microsatellite markers, ITS-RFLP (internal transcribed spacer-restriction fragment length polymorphism), RAPD (random amplified polymorphic DNA), and 2-DGE (two-dimensional gel electrophoresis) (e.g., Blok et al., 1997; Grenier et al., 2001; Plantard et al., 2008; Boucher et al., 2013). However, these methods focus only on specific sections of the genome, yield very few markers, or do not allow precise comparisons among populations. For example, Eves-van den Akker et al. (2015) recently published a metagenetic approach to genotype populations of the pale cyst nematode, *Globodera pallida*. While this method was shown powerful to rapidly evaluate genetic diversity and distribution of specific mitotypes, it is based on very few neutral markers, which prevent any assessment for selection. Genome scan is an interesting approach for the identification of genetic loci involved in adaptation to specific selection pressure. It was notably used in *G. pallida* to identify genomic regions associated with virulence on resistant potato cultivars (Eoche-Bosy et al., 2017).

Despite rapid advances in next-generation sequencing (NGS) technologies, sequencing a large number of individual genomes at high coverage in order to perform population genetic studies remains very expensive and may require important quantities of DNA from individuals. Elshire et al. (2011) developed a genotyping-by-sequencing (GBS) protocol to rapidly identify single-nucleotide polymorphisms (SNPs). The GBS technique and other restriction-site-associated DNA (RAD) sequencing methods produce large quantities of reads that do not cover the entire genome but have higher sequencing depth, thus reducing sequencing errors (Gautier et al., 2013; Anand et al., 2016). Loci generated by GBS can be present in both coding and non-coding regions and will be shared between all populations owing to the conservation of restriction sites (Cariou et al., 2013). Finally, GBS does not require any prior genomic information for the species being studied, which is an important consideration for SCN since there is no reference genome yet. The optimal gene coverage to reduce the amount of missing data depends on the choice of restriction enzyme (Fu et al., 2016). Fortunately, optimal gene coverage was already tested for the closely related species *Globodera rostochiensis* (Mimee et al., 2015).

The main objectives of this study were to (1) investigate the genetic relationships among SCN populations from United States
and Canada, (2) detect isolation by geographical distance (IBD) between SCN populations from United States and Canada, (3) detect genetic loci under selection associated to environmental and climatic parameters, and (4) identify the putative gene functions contributing to the adaptation of SCN populations to specific environmental conditions.

# MATERIALS AND METHODS

## Soybean Cyst Nematode Populations
Sampling and Genotyping-by-Sequencing

A total of 64 field populations of SCN, representative of the area currently infested in North America, were sampled or provided by collaborators from 11 US states (Delaware, Iowa, Illinois, Indiana, Kansas, Michigan, Minnesota, Missouri, North Dakota, Ohio, South Dakota) and one Canadian province (Ontario) (Figure 1). For DNA extraction, 40 cysts were randomly chosen from each population and pooled together. Eggs were extracted from each cyst and then washed twice in sterile filtered water. Total genomic DNA of each pool was extracted using the DNeasy Blood and Tissue Kit (Qiagen, Mississauga, ON, Canada) following the manufacturer's instructions. DNA extracts were quantified using Qubit fluorometric quantification (ThermoFisher Scientific, Burlington, ON, Canada) and normalized at 2 ng/µL prior to library preparation and sequencing. These steps were performed following standard protocols (Elshire et al., 2011; Poland et al., 2012) at the Institute for Integrative and Systems Biology (IBIS; Université Laval, Québec City, QC, Canada). Genotyping-by-sequencing was performed using the method described by MImmee et al. (2015) with a combination of two restriction enzymes (PstI/MspI) (New England Biolabs, Whitby, ON, Canada). After the restriction enzyme treatment, all samples (one composite per field population) were barcoded and multiplexed to obtain a single library for the 64 samples, which was sequenced on three Ion Proton chips (ThermoFisher Scientific) at IBIS.

## Single-Nucleotide Polymorphism Calling

The UNEAK pipeline (Lu et al., 2013), which is part of the TASSEL 3.0 bioinformatics analysis package (Bradbury et al., 2007), was used to process the raw reads, since no complete reference genome is yet available for *H. glycines*. This pipeline is designed to call SNPs *de novo*, without a reference genome with high stringency. Only sequences containing a single putative SNP (1-bp mismatch) per sequencing read were kept. Before analysis, the SNPs were further filtered with a minimum coverage (minCov) of 20 reads, a maximum coverage (maxCov) of 10,000 reads, and a minimum allele frequency (MAF) of 0.01.

## Population Genetics

Clustering of SCN populations using PCA was carried out using the *prcomp* function from the *stats* package in R software (R Core Team, 2017). The *poppr* package (v2.4.1) in R (Kamvar et al., 2014) was used to investigate the genetic relationship between populations and to build a phylogenetic tree (Provesti's distance, 10,000 bootstraps, neighbor-joining algorithm counting missing data as equivalent in the distance computation). Visualization of

sample coordinates and phylogenetic relationship analyses were carried out using the *phytools* 0.6–20 package in R (Revell, 2017).

Fixation index ($F_{ST}$) values (Wright, 1943) were calculated using the classical approach (Hartl and Clark, 2007) with the PoPoolation2 software (Kofler et al., 2011) to evaluate the genetic differentiation between each pair of populations. Furthermore, the effect of isolation by distance on population structure was tested using the correlation between the genetic distance ratio [$F_{ST}$/(1-$F_{ST}$)], as defined by Rousset (1997), and the geographic distance of population pairs in kilometers (km). The geographic distance between each sample location was calculated with the haversine formula using the geosphere package in R. To examine the significance of the relationship between the genetic distance ratio and the pairwise geographic distance distances, we performed a Mantel test (Spearman rank correlation) using the vegan package in R with 999 permutations. To evaluate the effect of time on population genetic differentiation, linear regression analyses were run in R software using pairwise $F_{ST}$ distances based on a point of origin selected on the basis of the first reports of SCN in North America. In our dataset, the closest sample to the oldest population was located in Clarkton, Missouri (MO1) (Hegge, 1957; Tylka and Marett, 2017). To access the genetic isolation from the MO1 location, a Spearman rank correlation ($r_s$) test was performed at a 0.95 confidence level.

Populations from three states (North Dakota, Minnesota, Illinois) were selected as case studies to evaluate the local (short-scale) genetic differentiation. These states were chosen because there were sufficient samples for comparison and because (i) North Dakota shows a recent introduction of SCN and a continuous northward dispersal of the nematode (Nelson and Bradley, 2003; Mathew et al., 2015); (ii) Minnesota has a longer history of SCN, with the nematode being first detected in 1978 and the infested area still expanding each year (Zheng et al., 2006); (iii) in Illinois SCN has been well established in every county for many years (Riggs, 2004; Tylka and Marett, 2017).

## Outlier Detection and Their Association
With Environmental Variables

Two geographic and four climatic covariables were investigated as possible factors explaining loci under selection: latitude (LAT), longitude (LONG), annual mean air temperature (BIO1), maximum air temperature of the warmest month (BIO5), annual precipitation (BIO12), and total precipitation of the warmest quarter (BIO18). All climatic variables were retrieved from the WorldClim global climate database, version 1.4 (Fick and Hijmans, 2017), corresponding to historic conditions (1960–1990). The spatial resolution used for the bioclimatic analysis was set at 30 s or 0.86 km².

To detect correlations between variations in population allele frequencies of SNPs and environmental factors, we used three different Bayesian methods (software programs): BayPass, version 2.1 (Gautier, 2015), BayeScan, version 2.1 (Foll and Gaggiotti, 2008), and BayeScEnv, version 1.1 (de Villemereuil and Gaggiotti, 2015). For each program, triplicate runs with different random seeds were performed with a pilot run of
<figure><img src="image_6.png" /><figcaption>**FIGURE 1** | Genetic relationship among North American populations of the soybean cyst nematode. (A) Three-dimensional principal component analyses (PCAs). (B) Phylogenetic tree plot based on Provesti's absolute genetic distance. Significant bootstrap values (> 50%) are indicated. Soybean cyst nematode populations from Ontario (Canada), North Dakota and Minnesota, and central US states are indicated in red, green and black, respectively.</figcaption></figure>

10,000 iterations to estimate starting parameters, a burn-in length of 50,000 iterations, and a minimum of 50,000 iterations, accounting for the small number of SNPs being investigated.

Post-run diagnostics were carried out using the *coda* package in R software in order to ensure sufficient iterations and normality of the Markov chain (Plummer et al., 2006).
The first software, BayPass, identifies genetic markers subject to selection by covariates such as phenotypic or environmental variables associated with the population of interest (Gautier, 2015). This application is based on the BAYENV model proposed by Coop et al. (2010) and Günther and Coop (2013), but with several modifications detailed in Gautier (2015), including the reprogramming of the Markov chain Monte Carlo (MCMC) algorithm. For BayPass, we used a pool-size file with the -d0yij option set at 800 and 20 pilot runs with a length of 10,000 iterations, a burn-in length of 50,000 iterations, and a chain length of 50,000 iterations. For each SNP, BayPass generates a Bayes factor (BF), quantifying evidence against the null hypothesis, and an empirical Bayesian *p*-value (eBPis), a metric measuring the difference between observed data and a simulated set of data (posterior distribution) (Kass and Raftery, 1995; Andraszewicz et al., 2015; Gautier, 2015). To be considered under selection, a SNP had to meet two criteria: a BF greater than 10 (BF > 10) and an eBPis lower than 0.05 (eBPis < 0.05).

The second software, BayeScan, uses a Bayesian likelihood method that assumes a Dirichlet distribution of allele frequencies between populations (Foll, 2012). This program estimates the probability that each locus is subject to selection by using a logistic regression on the two locus-population $F_{ST}$ coefficients. This Bayesian method uses a reversible-jump MCMC algorithm to calculate a posterior probability that each locus is under selection. The decision criterion to determine whether a locus is likely to be under a strong selection is the *q*-value (Foll, 2012), analog to a false discovery rate (FDR) *p*-value, that must be under 0.05. A second decision criterion was applied to further endorse the selected outliers: the ratio of posterior probabilities (PO). The PO threshold to affirm that a locus was under selection, in comparison with a neutral model, was set to 0.91, which corresponds to a strong relationship on the Jeffreys scale (Foll, 2012). We used BayeScan with the default parameters, but we set the minimum number of iterations to 50,000, the length of 20 pilot runs to 10,000 iterations, and the burn-in length to 50,000 iterations.

The third software, BayeScEnv, is similar to BayeScan and uses the $F_{ST}$ index to detect loci with a high level of differentiation in comparison with the entire genome. This program allows a normalization vector to be applied to the environmental data instead of only a binary combination, thus generating a lower number of false positives, according to the creators of the software (de Villemereuil and Gaggiotti, 2015). We used BayeScEnv with the default parameters, with the number of iterations set to 50,000 and 20 pilot runs with a length of 10,000 iterations. We used the reported *q*-value, which is related to the FDR, as our decision criterion, considering only SNPs with a *q*-value less than 0.05.

## Localization of Outlier Loci, Gene Function, and Genic Environment

The SNPs identified in short reads were first retrieved from a draft SCN genome available from SCNBase (https://www.scnbase.org/) by means of BLASTN with the default parameters, except for a smaller word size of 4, with the Blast2GO application

(Conesa et al., 2005). Many of the identified SNP-containing fragments matched multiple genes or genome locations (see Supplementary Tables 3–5). To assign a putative gene function to each SNP, we compared the aligned sequences to the National Center for Biotechnology Information (NCBI) protein database by means of BLASTX and BLASTP (Altschul et al., 1990) on a subset of sequences (nematodes, taxid: 6231) or to all of the NCBI non-redundant (nr) sequence database with an *E*-value significance cutoff of $10^{-5}$.

As our GBS sequencing covers approximately 0.8% of the genome (see Results section), we explored the genomic regions around outlier loci for genetic hitchhiking in whole-genome sequences from four populations of different origin and distinct virulence profile (ON1, ON34, IL4, and KS2). Sequencing libraries were prepared using the Nextera XT DNA Library Preparation Kit (Illumina, San Diego, CA) and sequenced on a MiSeq sequencer (Illumina) using the MiSeq Reagent Kit v3 (600-cycle). Reads were demultiplexed using the Sabre software program (https://github.com/najoshi/sabre) and processed with Trimmomatic (v0.32) (Bolger et al., 2014) to remove adapters and barcodes. Alignment on the draft reference genome (see above) was done using SAMtools (v0.1.19) (Li et al., 2009) and BWA (v0.7.12) (Li and Durbin, 2009). Only a window of ±50 kb around the 15 loci under selection was retained for the genic environment analysis. Variants were called with freebayes (v1.0.2) (Garrison and Marth, 2012) and snpeff (v4.2) (Cingolani et al., 2012). Distribution of allele frequencies in these four populations for each gene variant was compared with the allele frequency of the corresponding loci under selection in the same population in the genotyping-by-sequencing dataset. Only SNPs exhibiting a similar allele frequency and having a high impact on the predicted protein product were retained in our analysis.

## RESULTS

### Genotyping-by-Sequencing

A total of 192,576,709 short reads were obtained from the sequencing of the DNA from the 64 SCN populations, following digestion by the PstI/MspI restriction enzymes. The *H. glycines* effective genome length was estimated to be 96,752,286 bp, with an average of one PstI restriction site every 12,688 bases (E. Lord, private communication, 2017). On the basis of the sequence size used by the UNEAK (Universal Network-Enabled Analysis Kit) pipeline (first 64 bp of each read), the expected horizontal coverage was approximately 0.8% of the total SCN genome, and the vertical coverage at each locus, considering the number of reads obtained, was $400 \times$. Before filtering, the UNEAK pipeline identified 3,172 variants. After filtering for low coverage and selecting only SNP variants, two datasets were generated. The first dataset contained 245 SNPs without missing data (loci sequenced in all populations), and the second contained 804 SNPs with missing data (Sequences and coverage in Supplementary Table 7).

### Population Structure

#### Principal Component Analysis

A principal component analysis (PCA) was performed using the dataset without missing data, which contained 245 SNPs. The
PCA plot of the 64 North American SCN populations revealed a geographically ordered pattern (Figure 1A). Overall, the Ontario populations showed the greatest dispersion, indicating they are more genetically differentiated, while the central regions of the US showed more clustered populations. Also, the group containing populations from Minnesota and North Dakota (in green in Figure 1A) was clearly different from the Ontario populations by the first axis, which explained 21% of the total variation. All populations from the central states, although they originated from many states and covered a much wider area, were less diverse and more clustered together.

## Phylogenetic Tree

To better understand the population structure of each population, we used the 804-SNP dataset to conduct a phylogenetic analysis based on Provesti's distance. The inferred neighbor-joining phylogenetic tree (Figure 1B) showed that most of the Ontario samples (in red in Figure 1B) were different from the rest of the North American populations. Some samples from Ontario were found to be different from each others and more similar to those of central states. Such a pattern can also be observed in the PCA analysis (Figure 1A).

## Genetic Differentiation

Overall, we observed an average value for Wright's fixation index ($F_{ST}$) of $0.15 \pm 0.07$. However, pairwise genetic distance was highly variable, with a minimum $F_{ST}$ of 0.005 and a maximum $F_{ST}$ of 0.53 (Supplementary Table 1). The highest $F_{ST}$ value was found between two of the most remote populations, ND6 from North Dakota and ON41 from Ontario, which were separately by approximately 1,720 km. On the other hand, the lowest $F_{ST}$ value was observed between ON3 and ON5, which were separated by only 17.9 km. In general, higher $F_{ST}$ values were found in distant populations, especially between populations from Minnesota (MN4 and MN5), North Dakota (ND3 and ND6), South Dakota (SD2), and Ontario (ON6, ON13, ON25, ON41, and ON46).

## Isolation by Distance

The analysis of the relationship between genetic distance [$F_{ST}$/(1-$F_{ST}$)] and geographic distance (km) for all pairs of populations revealed a significant ($p = 0.012$) but weak effect of isolation by distance, with a Mantel correlation ($r$) of only 0.135. However, when analyzing the effect of isolation by distance at the local scale, a different pattern was found. With respect to US states with recent SCN history, a strong and significant positive correlation was observed between the genetic and geographic distances in North Dakota ($r = 0.430, p = 0.013$) and Minnesota ($r = 0.730, p = 0.006$). In contrast, no correlation was found among the populations from Illinois ($r = -0.170, p = 0.767$).

To further explore the contribution of isolation by distance to the genetic differentiation of United States and Canada populations of SCN, we compared the genetic distance of each population from a theoretical ancestor population (MO1) using the 245-SNP dataset. This analysis included 38 independent populations that were separated by 273 to 1,367 km from MO1 (Figure 2A). The genetic-distance ratio to the “oldest population” varied from 0.06 for IL1 (386 km apart) to 0.52 for

ND6 (1,327 km apart) (Supplementary Table 2) and an overall moderate ($r_s = 0.535$) but significant ($p = 0.001$) relationship was observed between genetic distance and geographic distance (Figure 2B).

## Signature of Genetic Adaptations

### Genome Scans

Using three Bayesian inference genome scan approaches, we identified 71 SNPs (out of 804) that were under selection. Among the different computational methods, 48 loci were identified by BayeScan, 23 loci were identified by BayPass, and 16 loci were identified by BayeScEnv (Supplementary Tables 3–5). Out of those 71 SNPs, only one was highlighted by the three methods, while 15 outliers were inferred by at least two of the methods (Figure 3). The geographic variables (LAT and LONG) were associated with seven and eight outlier loci, respectively, and the climatic variables (BIO1, BIO5, BIO12, and BIO18) were associated with nine, six, eight, and seven SNPs, respectively (Table 1). Three SNPs were strictly associated with

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>**FIGURE 3** | Number of outlier loci identified by three different Bayesian genome scan methods in 64 populations of the soybean cyst nematode from North America.</figcaption></figure>

only one environmental variable: TP188858 was associated with maximum temperature, while TP252086 and TP364091 appeared to be linked to annual temperature. The allele frequencies analysis at these 15 loci revealed SNPs that were specific for a given region, such as TP364091, TP181514, and TP333913 for Minnesota and North Dakota and TP376687, TP146577, and TP226227 for Ontario (**Figure 4**).

## Localization, Effect, and Genic Environment of Outlier Loci

Of the 15 SNPs under selection identified by at least two Bayesian approaches, 10 SNPs were located in predicted genes (seven in exons and three in introns), and five SNPs were located in intergenic regions (**Table 2**). Two SNPs (TP380188 and TP364091) correspond to non-synonymous mutations that induce a change in the amino acid sequence of the resulting protein. The annotations of the genes putatively impacted by these modifications are SMAD-4 (mothers against decapentaplegic homolog 4) for TP380188 and DOP-3 (dopamine receptor 3) for TP364091. The sequence corresponding to TP364091 was retrieved in four different genes (g14639.t1, g14642.t1, g14644.t1, and g14654.t1) in the draft reference genome and introduced two different mutations (Supplementary Tables 3–5). All the other SNPs (13/15) were silent mutations (synonymous in exons or located in introns or intergenic regions) and thus probably not the cause of their selection. Among the 42 remaining SNPs that were under selection but identified by only one pipeline, 12 correspond to a

non-synonymous modification impacting the predicted protein sequence. These SNPs include the gene UFM1 (Ubiquitin-fold modifier 1), the gene PLA2 (85/88 kDa calcium-independent phospholipase A2), and a pyruvate kinase (Supplementary Tables 3–5).

The whole-genome resequencing of four SCN populations for the exploration of the genic environment of outlier loci highlighted 257 genes containing high-impact variants in a ± 50-kb window around loci under selection. Annotation was available for 148 of these genes, corresponding to 132 different gene functions (Supplementary Table 6). Of these genes, 25 exhibited allele frequencies similar to their neighboring outlier loci (**Table 3**).

## DISCUSSION

Since its introduction into North America in 1954, the soybean cyst nematode has spread at a regular pace to most of the soybean-producing areas in the US and reached Canada in 1987. Its northern limit of establishment remains to be determined, as new populations are found each year. Soybean acreage is rapidly increasing in northern latitudes; for example, farmers in the province of Manitoba, Canada, sowed 2.3 million acres with soybean in 2017, a 40% increase in 1 year in a province where soybean was not grown a decade ago (Statistics Canada, 2017). This spectacular change in land use results from a coordinated breeding program for short-season and drought-tolerant soybean cultivars (Tardivel et al., 2014). Although soybean is nowadays very profitable in Canada (Manitoba, Northern Ontario, Québec), managing a pest such as SCN could be challenging, as no resistant cultivars are currently adapted to these regions. Simulations based on thermal development have shown that SCN could theoretically survive in these new soybean growing areas (Gendron St-Marseille, 2013), but actual establishment is a different matter, as initial observations in the province of Québec indicated that even though the nematode is detected in many regions, populations weakly reproduced (Mimee et al., 2016). This phenomenon probably reflects SCN's very recent history in this part of the world and the poor fitness of the introduced populations but does not exclude possible adaptations of the nematode in the future. To better understand the evolution of SCN population genetics in North America, we compared 64 populations originating from 11 US states and the province of Ontario in Canada.

The dynamics of cyst nematode populations are intrinsically complex. Several thousand cysts from different females can form on a single plant, and each cyst contains hundreds of half-sibling individuals as a result of SCN's polyandrous mode of reproduction. Thus, the number of unique genotypes occurring in a field is massive. Isolating individual nematodes, even hundreds of them, to explore the genetic relationship between populations would be an arduous task and the approach would still be biased owing to limited sub-sampling. For these reasons, we opted for a Pool-Seq (sequencing of pooled DNA samples) approach. Sequencing DNA from pooled samples for each population also has the advantage of keeping the number of
<table><caption>**TABLE 1** | Association of outlier loci [single-nucleotide polymorphisms (SNPs)] identified by at least two Bayesian genome scan methods in 64 populations of the soybean cyst nematode from North America with geographic or climatic co-variables.</caption><thead><tr><th>SNP ID</th><th>Method</th><th>Latitude</th><th>Longitude</th><th>Temperature (annual)</th><th>Temperature (maximum)</th><th>Precipitation (annual)</th><th>Precipitation (warmest quarter)</th></tr></thead><tbody><tr><td rowspan="3">TP66981</td><td>BayeScan</td><td></td><td>0.998</td><td>0.995</td><td></td><td></td><td></td></tr><tr><td>BayeScEnv</td><td></td><td>0.017</td><td>0.042</td><td></td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP77825</td><td>BayeScan</td><td></td><td></td><td></td><td>0.999</td><td>0.941</td><td></td></tr><tr><td>BayeScEnv</td><td></td><td></td><td></td><td>0.029</td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP146577</td><td>BayeScan</td><td>0.952</td><td>1.000</td><td></td><td>1.000</td><td>0.997</td><td>0.997</td></tr><tr><td>BayeScEnv</td><td></td><td>0.047</td><td></td><td>0.038</td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP173057</td><td>BayeScan</td><td>0.988</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>BayeScEnv</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td>13.86</td></tr><tr><td rowspan="3">TP181514</td><td>BayeScan</td><td>0.996</td><td>0.940</td><td>0.992</td><td>0.961</td><td></td><td></td></tr><tr><td>BayeScEnv</td><td></td><td>0.050</td><td>0.009</td><td></td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP188858</td><td>BayeScan</td><td></td><td></td><td></td><td>0.999</td><td></td><td></td></tr><tr><td>BayeScEnv</td><td></td><td></td><td></td><td>0.050</td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP226227</td><td>BayeScan</td><td></td><td>1.000</td><td></td><td>1.000</td><td>0.956</td><td>0.999</td></tr><tr><td>BayeScEnv</td><td></td><td></td><td></td><td>0.044</td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP252086</td><td>BayeScan</td><td></td><td></td><td>0.963</td><td></td><td></td><td></td></tr><tr><td>BayeScEnv</td><td></td><td></td><td>0.050</td><td></td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP324827</td><td>BayeScan</td><td>0.999</td><td>0.999</td><td>0.999</td><td></td><td>1.000</td><td>0.999</td></tr><tr><td>BayeScEnv</td><td>0.043</td><td>0.023</td><td>0.006</td><td></td><td>0.017</td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP333577</td><td>BayeScan</td><td>0.984</td><td>0.978</td><td></td><td></td><td>0.999</td><td>1.000</td></tr><tr><td>BayeScEnv</td><td></td><td>0.029</td><td>0.015</td><td></td><td></td><td>0.040</td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td>14.32</td></tr><tr><td rowspan="3">TP333913</td><td>BayeScan</td><td>0.931</td><td>0.978</td><td>0.976</td><td></td><td>0.930</td><td></td></tr><tr><td>BayeScEnv</td><td></td><td>0.029</td><td>0.015</td><td></td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP364091</td><td>BayeScan</td><td></td><td></td><td>0.976</td><td></td><td></td><td></td></tr><tr><td>BayeScEnv</td><td></td><td></td><td>0.020</td><td></td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP376687</td><td>BayeScan</td><td>0.993</td><td>0.999</td><td></td><td>1.000</td><td>0.999</td><td>0.999</td></tr><tr><td>BayeScEnv</td><td></td><td></td><td></td><td>0.038</td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP380020</td><td>BayeScan</td><td></td><td></td><td>0.998</td><td></td><td>1.000</td><td></td></tr><tr><td>BayeScEnv</td><td></td><td></td><td>0.030</td><td></td><td>0.027</td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">TP380188</td><td>BayeScan</td><td></td><td></td><td>0.968</td><td></td><td></td><td></td></tr><tr><td>BayeScEnv</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>BayPass</td><td></td><td></td><td></td><td></td><td></td><td>13.97</td></tr></tbody></table>

redundant sequences low (Schlötterer et al., 2014). However, this approach does not allow the assignment of sequences to individuals and prevents some population genetic analyses such

as linkage-disequilibrium estimation (Lynch et al., 2014; Anand et al., 2016). Nevertheless, a number of recent studies have overcome these limitations and have proposed new methods to
<figure><img src="image_9.png" /><figcaption>**FIGURE 4** | Allele frequencies of 15 outlier loci identified in North American populations of the soybean cyst nematode and associations of the loci with geographic or climatic co-variables. Black squares represent a significant association.</figcaption></figure>

analyze genetic variations among populations based on Pool-Seq data (Van Tassell et al., 2008; Gautier et al., 2013; Navon et al., 2013; M, 2015; Anand et al., 2016).

Both the PCA and phylogenetic analyses showed strong clustering of populations based on their geographic origin, which supports the “nearest-neighbor theory” movement as in the stepping-stone model (Kimura and Weiss, 1964; Hutchison and Templeton, 1999). We observed a geographic separation dividing SCN populations into two clusters, a northeastern one (Ontario) and a northwestern one (Minnesota, North Dakota) that both share genetic similarities with central populations. In northeast and northwest regions, SCN populations established only recently, and new detections continue to occur in the northern areas. In the central states, SCN populations have been present for a longer time and resistant soybean cultivars are routinely used. We believed that these two conditions contribute to explain the homogeneity of the central populations. On one hand, the continuous exchange of genetic material would have led to the homogenization of alleles in sympatric populations. On the

other hand, the massive use of nematode resistant cultivars that are all derived from a single source of resistance (PI 88788) would have selected for the virulent SCN genotypes, thereby limiting the propagation of some alleles. It is thus surprising to observe a greater genetic differentiation in SCN populations from Ontario since the main assumption is that they originate from central states. These populations should contain less diversity, owing to the founder effect, which would normally result in reduced potential for differentiation. However, the selective pressures mentioned above are not yet present in Ontario, allowing the development of most genotypes. This finding concurs with the results of Faghihi et al. (2010), who observed a much greater diversity in SCN population phenotypes (HG types) in Ontario than in Tennessee, Indiana, or Illinois. A similar pattern was obtained in Minnesota by Zheng et al. (2006) and in China by Liu et al. (1997) when comparing populations from the north to populations from central regions. However, all these authors suggested that climatic conditions or local environmental factors may also play a role.
<table><caption>**TABLE 2** | Identification and effect of outlier loci [single-nucleotide polymorphisms (SNPs)] identified by at least two Bayesian genome scan methods in 64 populations of the soybean cyst nematode from North America.</caption><thead><tr><th>ID</th><th>Location</th><th>S or N-S<sup>1</sup></th><th>SNP</th><th>a.a. subst.<sup>2</sup></th><th>E-value</th><th>Best annotation</th></tr></thead><tbody><tr><td>TP380188</td><td>Exon</td><td>N-S</td><td>C/A</td><td>Gln/Lys</td><td>2.0e-20</td><td>SMAD-4 [Strongyloides ratti]</td></tr><tr><td>TP364091</td><td>Exon</td><td>N-S</td><td>G/A</td><td>Asp/Asn</td><td>2.6e-104</td><td>Dopamine receptor 3 [Toxocara canis]</td></tr><tr><td>TP376687</td><td>Exon</td><td>S</td><td>T/G</td><td>Gly/Gly</td><td>4.0e-11</td><td>Putative esophageal gland cell protein Hgg-20 [Heterodera glycines]</td></tr><tr><td>TP173057</td><td>Exon</td><td>S</td><td>C/T</td><td>Ile/Ile</td><td>2.0e-11</td><td>Leucine-rich repeat-containing domain protein, partial [Necator americanus]</td></tr><tr><td>TP333577</td><td>Exon</td><td>S</td><td>C/T</td><td>Leu/Leu</td><td>1.0e-110</td><td>Cadherin-4 [Caenorhabditis elegans]</td></tr><tr><td>TP181514</td><td>Exon</td><td>S</td><td>T/C</td><td>Gly/Gly</td><td>1.6e-100</td><td>LET-19 [Caenorhabditis brenneri]</td></tr><tr><td>TP333913</td><td>Exon</td><td>S</td><td>T/C</td><td>Gly/Gly</td><td>1.6e-7</td><td>EGL-17 [Brugia malayi]</td></tr><tr><td>TP252086</td><td>Intron</td><td>–</td><td>G/A</td><td>–</td><td>2.0e-64</td><td>Esophageal gland-localized secretory protein 11 [Heterodera glycines]</td></tr><tr><td>TP324827</td><td>Intron</td><td>–</td><td>C/T</td><td>–</td><td>3.0e-32</td><td>Reverse transcriptase [Ancylostoma ceylanicum]</td></tr><tr><td>TP146577</td><td>Intron</td><td>–</td><td>T/A</td><td>–</td><td>2.0e-11</td><td>Hypothetical protein Y032_0209g2110 [Ancylostoma ceylanicum]</td></tr><tr><td>TP77825</td><td>Intergenic</td><td>–</td><td>T/A</td><td>–</td><td>2.0e-57</td><td>Nearest gene: hypothetical protein CRE_19487 [Caenorhabditis remanei]</td></tr><tr><td>TP226227</td><td>Intergenic</td><td>–</td><td>T/G</td><td>–</td><td>3.0e-05</td><td>Nearest gene: ATP-dependent DNA helicase PIF1, partial [Toxocara canis]</td></tr><tr><td>TP66981</td><td>Intergenic</td><td>–</td><td>T/G</td><td>–</td><td>–</td><td>No significant hit</td></tr><tr><td>TP380020</td><td>Intergenic</td><td>–</td><td>G/T</td><td>–</td><td>–</td><td>No significant hit</td></tr><tr><td>TP188858</td><td>Intergenic</td><td>–</td><td>G/A</td><td>–</td><td>–</td><td>No significant hit</td></tr></tbody></table>

According to Hartl and Clark (2007), the overall average $F_{ST}$ value of 0.15 that we obtained by comparing all pairs of SCN populations from North America corresponds to a high level of differentiation. However, when we examined specific pairs of populations, we found that the populations were very similar at the local scale. Also, the differentiation of recent populations (Ontario, Minnesota, North Dakota) from those of the central states was only moderate, whereas comparing these recent populations together generated very high $F_{ST}$ values (up to 0.53). This corroborates results from the PCA analysis and clearly suggests that populations diverged following two distinct routes, one in the northwest toward North Dakota and one in the northeast toward Ontario. This hypothesis is further supported by the isolation-by-distance analysis. Overall, we observed a significant effect of distance on genetic differentiation at the continental scale. The effect of distance was even more noticeable in newly colonized areas like North Dakota and Minnesota or when each population was compared to a common ancestor. On the other hand, the absence of any effect of distance in Illinois and low $F_{ST}$ in the central states clearly suggests that dispersal is frequent enough to weaken signs of isolation by distance after a short period of time. This continuous gene flow combined with the selection of alleles due to the general use of the same resistant soybean lines are thought to cause a rapid homogenization of SCN populations (Mitchum et al., 2007; Niblack et al., 2008; Zheng and Chen, 2011). If SCN is already adapted to environmental conditions prevailing in northern areas, this homogenization process would accelerate the development of virulent genotypes in these regions (in comparison with the central states) and would be independent from the use of resistant cultivars. This concurs with the finding of Faghihi et al. (2010) where 15% of SCN populations in Ontario were able to multiply on PI 548402 and 6% were able to multiply

on PI 90763, even though these sources of resistance were not present in commercial cultivars in Ontario.

Although 71 SNPs were found to be under selection among SCN populations, only one SNP was identified by the three Bayesian approaches tested. This SNP (TP333577) was strongly associated with precipitation; but the geographic distribution of allele frequencies (Figure 4) at that specific locus did not match with the general pattern of genetic differentiation, instead showing a very local effect (in populations from Minnesota but also in some populations from Iowa, Illinois, North Dakota, and Ohio). Thus, this selection for wetness was either lost in favor of a more stringent pressure or too costly to be maintained in the other regions. Overall, considering the broad range of environmental conditions in which these populations develop, finding only a few genes under selection suggests that SCN possesses an intrinsic capacity to evolve within a large range of temperature and wetness gradients.

Most of the 15 outliers identified by at least two pipelines correlated with temperature, reinforcing that temperature is a key factor that acts on the genetic diversity and natural selection of nematode populations. Nematodes are poikilothermic organisms and depend greatly on temperature to complete their life cycle. Thus, any adaptation that facilitates the development and reproduction of nematodes at lower temperatures would be advantageous. On the basis of the distribution of allele frequencies, however, this adaptation does not seem crucial for survival, as the northern populations did not all exhibit similar patterns. This trait is probably not yet fixed in these populations. Nevertheless, six SNPs under selection and associated with temperature were specific to populations from Ontario (TP376687, TP146577, and TP226227) or only retrieved from Minnesota and North Dakota (TP364091, TP181514, and TP333913), a finding that clearly suggests local adaptations.

<footnote n="1">¹Synonymous (S) or non-synonymous (N-S) mutations in the final protein product.</footnote>

<footnote n="2">²a.a., amino acid substitution.</footnote>
<table><caption>**TABLE 3** | Putative hitchhiking gene variants with high-impact neighboring loci under selection and exhibiting a similar allele frequency in *Heterodera glycines*.</caption><thead><tr><th>ID</th><th>Gene</th><th>Best annotation</th><th>E-value<sup>1</sup></th><th>Type of variants<sup>2</sup></th></tr></thead><tbody><tr><td>TP146577</td><td>g15502.t1</td><td>G-protein-coupled receptor, partial [Pristionchus pacificus]</td><td>9.10e-57</td><td>complex, insertion</td></tr><tr><td>TP173057</td><td>g18365.t1</td><td>85 kda calcium-independent phospholipase a2 [Ascaris suum]</td><td>8.50e-46</td><td>SNP, insertion</td></tr><tr><td>TP173057</td><td>g18358.t2</td><td>F-box only protein 30 [Ascaris suum]</td><td>1.10e-33</td><td>4 complexes</td></tr><tr><td>TP173057</td><td>g18327.t1</td><td>Protein CBR-CLEC-223 [Caenorhabditis briggsae]</td><td>5.70e-57</td><td>SNP</td></tr><tr><td>TP173057</td><td>g18328.t1</td><td>Protein lap1 [Toxocara canis]</td><td>1.80e-104</td><td>6 SNP, 3 insertions, 2 complexes</td></tr><tr><td>TP173057</td><td>g18332.t1</td><td>Protein scribble-like protein, partial [Stegodyphus mimosarum]</td><td>1.00e-33</td><td>2 deletions, MNP</td></tr><tr><td>TP181514</td><td>g3633.t1</td><td>Mediator of RNA polymerase II transcription subunit 13-like protein [Ascaris suum]</td><td>1.20e-133</td><td>Deletion</td></tr><tr><td>TP226227</td><td>g1982.t1</td><td>85/88 kDa calcium-independent phospholipase A2 [Toxocara canis]</td><td>8.00e-29</td><td>SNP</td></tr><tr><td>TP226227</td><td>g14545.t1</td><td>Glycoside hydrolase domain containing protein [Haemonchus contortus]</td><td>1.10e-162</td><td>deletion, deletion</td></tr><tr><td>TP226227</td><td>g12205.t1</td><td>Kazal domain and organic anion transporter polypeptide (OATP) family and MFS [Strongyloides ratti]</td><td>3.00e-36</td><td>deletion, 3 complexes</td></tr><tr><td>TP226227</td><td>g1979.t1</td><td>Serine/threonine-protein kinase haspin [Haemonchus contortus]</td><td>1.40e-62</td><td>insertion</td></tr><tr><td>TP252086</td><td>g18812.t1</td><td>Bestrophin/UPF0187 family-containing protein [Strongyloides ratti]</td><td>2.30e-44</td><td>deletion</td></tr><tr><td>TP252086</td><td>g18825.t1</td><td>mtN3/saliva family protein [Oesophagostomum dentatum]</td><td>1.30e-55</td><td>deletion, SNP, complex</td></tr><tr><td>TP333577</td><td>g16696.t2</td><td>AFG3-like protein 2 [Toxocara canis]</td><td>1.30e-286</td><td>complex, 2 insertions</td></tr><tr><td>TP333913</td><td>g14633.t1</td><td>Cytohesin-1 [Trichinella zimbabwensis]</td><td>1.30e-118</td><td>SNP</td></tr><tr><td>TP333913</td><td>g14626.t1</td><td>Fibroblast growth factor 17 [Toxocara canis]</td><td></td><td>deletion</td></tr><tr><td>TP333913</td><td>g14617.t1</td><td>Immunoglobulin domain containing protein [Haemonchus contortus]</td><td>9.60e-42</td><td>deletion</td></tr><tr><td>TP333913</td><td>g14614.t1</td><td>Sodium-dependent acetylcholine transporter [Toxocara canis]</td><td>1.50e-215</td><td>SNP</td></tr><tr><td>TP364091</td><td>g14633.t1</td><td>Cytohesin-1 [Trichinella zimbabwensis]</td><td>1.30e-118</td><td>SNP</td></tr><tr><td>TP364091</td><td>g14656.t1</td><td>Major facilitator superfamily MFS-1 domain containing protein, partial [Haemonchus contortus]</td><td>2.00e-106</td><td>deletion</td></tr><tr><td>TP380020</td><td>g7291.t1</td><td>Carboxypeptidase [Eupolyphaga sinensis]</td><td>1.00e-60</td><td>SNP</td></tr><tr><td>TP380020</td><td>g8313.t1</td><td>Cell division cycle related [Caenorhabditis elegans]</td><td>1.10e-254</td><td>deletion, SNP</td></tr><tr><td>TP380020</td><td>g8311.t1</td><td>GPI ethanolamine phosphate transferase 3, partial [Trichinella papuae]</td><td>1.90e-125</td><td>deletion</td></tr><tr><td>TP380020</td><td>g8314.t3</td><td>Homeobox protein ceh-6 [Toxocara canis]</td><td>4.60e-52</td><td>insertion, 2 deletions</td></tr><tr><td>TP380020</td><td>g8309.t1</td><td>Immunoglobulin I-set domain containing protein [Ascaris suum]</td><td>5.50e-31</td><td>2 deletions</td></tr><tr><td>TP380020</td><td>g8305.t1</td><td>RecQ-mediated genome instability protein 1 [Toxocara canis]</td><td>4.60e-28</td><td>SNP</td></tr><tr><td>TP380188</td><td>g2288.t1</td><td>Cleavage polyadenylation specificity factor domain containing protein [Haemonchus contortus]</td><td>3.90e-209</td><td>insertion</td></tr><tr><td>TP380188</td><td>g2287.t1</td><td>Cytochrome b-c1 complex subunit 2, mitochondrial [Strongyloides ratti]</td><td>1.20e-63</td><td>deletion, SNP</td></tr></tbody></table>

The SNP corresponding to TP364091 was retrieved in four out of six populations from Minnesota and in one population from North Dakota. That SNP, which is a non-synonymous mutation in the gene DOP-3 (dopamine receptor 3) that replaces an aspartic acid with an asparagine, was strongly associated with the annual mean temperature and annual precipitation. This modification was retrieved in two distinct gene sequences (g14642.t1 and g14654.t1), and a synonymous mutation (Pro/Pro) was retrieved in two other sequences (g14639.t1 and g14644.t1) but in a different amino acid, indicating a reading frame shift and a putative upstream indel. These four genes were all annotated as DOP-3 and were probably the result of gene duplication. Gene duplication has frequently been hypothesized to play an important role in adaptation to the environment (reviewed in Kondrashov, 2012).

In *Caenorhabditis elegans*, the *DOP-3* protein is well known to modulate chemosensory functions, such as mating and foraging, and to be involved in locomotion (Chase et al., 2004; Wood and Ferkey, 2016). The other SNP that introduces a non-synonymous mutation and was highlighted by two genome scan methods, namely, TP380188, was located in the gene SMAD-4 (mothers against decapentaplegic homolog 4) and was retrieved in only a few populations from Minnesota, North Dakota, and Ontario. The SMAD-4 protein is a regulator of different cellular processes, including cell differentiation, apoptosis, migration, and proliferation (Nikolic et al., 2011). Mutations of SMAD-2, SMAD-3, SMAD-4, or SMAD-6 were previously shown to result in a 30% reduction in body size in *C. elegans* (Savage-Dunn et al., 2003; Watanabe et al., 2007). Also of interest is the locus TP122360, identified only by BayPass, that was located

<footnote n="1">¹Expected value for annotations of each putative hitchhiking gene.</footnote>

<footnote n="2">²Variant in bold indicates a stop codon gain. More details on types of variants and position in the gene are provided in Supplementary Table 6. SNP, single-nucleotide polymorphism; MNP, multiple-nucleotide polymorphism.</footnote>
in the UFM1 (ubiquitin-fold modifier 1) genes (g7250.t1 and g7228.t1). A deletion in this gene in *C. elegans*, although reducing reproduction rate and life span, increased the survival of this nematode under oxidative or heat stress (Hertel et al., 2013).

For the majority of the loci, their selection probably results from genetic hitchhiking rather than a direct contribution to adaptation. To explore that possibility, we analyzed the genic environment in a ± 50 kb window around each SNP in four SCN populations. Genes containing high-impact genetic variations and exhibiting the same allele frequencies as the associated SNP under selection were of particular interest. Two of these genes (g18327.t1 and g15502.t1) coded for G-protein-coupled receptors (GPCR). These proteins are members of a large and very diverse multigene family with hundreds of occurrences in the *C. elegans* genome (Bockaert and Pin, 1999; Bargmann, 2006). The GCPR are crucial in sensing the local environment and were shown to evolve following alterations in habitat or foraging behavior (Nei et al., 2008). Four genes coding for proteins implicated in metabolite transport into the cell were also identified (g12205.t1, g14656.t1, g18812.t1, and g18825.t1), as was one gene involved in the regulation of those proteins (g18365.t1). All these proteins are required to maintain homeostasis in the cell and to respond to local environmental changes. In a similar study using the fungus *Fagus sylvatica*, Pluess et al. (2016) found that a version of a potassium transporter was associated with lower precipitation and could contribute to the regulation of growth under dry conditions. Structural changes in the body of an organism as a result of microevolution can also confer a significant advantage in terms of resisting more adverse environmental conditions, such as drought or high temperatures (Hazel and Williams, 1990). In our study, mutations were observed in two genes involved in cell-membrane and cell-wall stability (g18328.t1 and g8311.t1). Lastly, modifications in genes involved in the regulation of transcription (g2288.t1 and g3633.t1) and in the maturation of proteins (g14633.t1 and g18358.t2) have the potential to radically change the proteome of the adapted organism.

Although some SCN genes were found to be under selection and local adaptation was found to be underway in this study, our results also indicate that there is no critical adaptive event required for SCN establishment in northern latitudes.

Consequently, all populations should theoretically survive and multiply at high latitudes. The risk is thus real for new soybean areas where cultivars resistant to SCN are not available at this time. Of course, the method we used does not explore the entire genome, and a pan-genomic study of these populations could reveal other loci under selection in an evolutionary process.

## DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

## AUTHOR CONTRIBUTIONS

BM, A-FG, and JB conceived and designed the experiments. A-FG, EL, and P-YV performed the experiments and analyzed the data. A-FG, EL, and BM wrote the paper. All the authors revised the manuscript.

## ACKNOWLEDGMENTS

The authors thank Drs. Thomas Baum and Andrew Severin for graciously making SCN genome sequences available to us before public release on scnbase.org. The authors also thank Tom Welacky (AAFC) and Terry Niblack (OSU) for providing SCN populations. The authors wish to acknowledge the contribution of Génome Québec, Genome Canada, the Government of Canada, the Ministère de l'Économie, Science et Innovation du Québec, Semences Prograin Inc., Syngenta Canada Inc., Sevita Genetics, La Coop Fédérée, Grain Farmers of Ontario, Saskatchewan Pulse Growers, Manitoba Pulse and Soybean Growers, the Canadian Field Crop Research Alliance and Les Producteurs de grains du Québec.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00987/full#supplementary-material

## REFERENCES

* Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic local alignment searchtool. *J. Mol. Biol.* 215, 403–410. doi: 10.1016/S0022-2836(05)80360-2

* Anand, S., Mangano, E., Barizzone, N., Bordini, R., Sorosina, M., Clarelli, F., et al. (2016). Next generation sequencing of pooled samples: guideline for variants’ filtering. *Sci. Rep.* 6:33735. doi: 10.1038/srep33735

* Anderson, T. R., Welacky, T. W., Olechowski, H. T., Ablett, G., and Ebsary, B. A. (1988). First report of *Heterodera glycines* on soybeans in Ontario, Canada. *Plant Dis.* 72:453. doi: 10.1094/PD-72-0453C

* Andraszewicz, S., Scheibehenne, B., Rieskamp, J., Grasman, R., Verhagen, J., and Wagenmakers, E.-J. (2015). An introduction to Bayesian hypothesis testing for management research. *J. Manag.* 41, 521–543. doi: 10.1177/0149206314560412

* Bargmann, C. I. (2006). Comparative chemosensation from receptors to ecology. *Nature* 444, 295–301. doi: 10.1038/nature05402

* Blok, V. C., Phillips, M. S., and Harrower, B. E. (1997). Comparison of British populations of potato cyst nematodes with populations from continental Europe and South America using RAPDs. *Genome* 40, 286–293. doi: 10.1139/g97-040

* Bockaert, J., and Pin, J. P. (1999). Molecular tinkering of G protein-coupled receptors: an evolutionary success. *EMBO J.* 18, 1723–1729. doi: 10.1093/emboj/18.7.1723

* Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for illumina sequence data. *Bioinformatics* 30, 2114–2120. doi: 10.1093/bioinformatics/btu170

* Boucher, A. C., Mimee, B., Montarry, J., Bardou-Valette, S., Bélair, G., Moffett, P., et al. (2013). Genetic diversity of the golden potato cyst nematode *Globodera rostochiensis* and determination of the origin of populations in Quebec, Canada. *Mol. Phylogenet. Evol.* 69, 75–82. doi: 10.1016/j.ympev.2013.05.020

* Bradbury, P. J., Zhang, Z., Kroon, D. E., Casstevens, T. M., Ramdoss, Y., and Buckler, E. S. (2007). TASSEL: software for association mapping
of complex traits in diverse samples. *Bioinformatics* 23, 2633–2635. doi: 10.1093/bioinformatics/btm308

Cariou, M., Duret, L., and Charlat, S. (2013). Is RAD-seq suitable for phylogenetic inference? An *in silico* assessment and optimization. *Ecol. Evol.* 3, 846–852. doi: 10.1002/ece3.512

Chase, D. L., Pepper, J. S., and Koelle, M. R. (2004). Mechanism of extrasynaptic dopamine signaling in *Caenorhabditis elegans*. *Nat. Neurosci.* 7, 1096–1103. doi: 10.1038/nn1316

Cingolani, P., Platts, A., Wang, L. L., Coon, M., Nguyen, T., Wang, L., et al. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SNPEff: SNPs in the genome of *Drosophila melanogaster* strain w¹¹¹⁸; *iso-2*; *iso-3*. *Fly* 6, 80–92. doi: 10.4161/fly.19695

Colgrove, A. L., and Niblack, T. L. (2008). Correlation of female indices from virulence assays on inbred lines and field populations of *Heterodera glycines*. *J. Nematol.* 40, 39–45.

Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. *Bioinformatics* 21, 3674–3676. doi: 10.1093/bioinformatics/bti610

Coop, G., Witonsky, D., Di Rienzo, A., and Pritchard, J. K. (2010). Using environmental correlations to identify loci underlying local adaptation. *Genetics* 185, 1411–1423. doi: 10.1534/genetics.110.114819

Davis, E. L., and Tylka, G. L. (2000). Soybean cyst nematode disease. *Plant Health Instr.* doi: 10.1094/PHI-I-2000-0725-01

de Villemereuil, P., and Gaggiotti, O. E. (2015). A new FST-based method to uncover local adaptation using environmental variables. *Meth. Ecol. Evol.* 6, 1248–1258. doi: 10.1111/2041-210X.12418

Early, R., and Sax, D. F. (2014). Climatic niche shifts between species' native and naturalized ranges raise concern for ecological forecasts during invasions and climate change. *Glob. Ecol. Biogeogr.* 23, 1356–1365. doi: 10.1111/geb.12208

Eizaguirre, C., and Baltazar-Soares, M. (2014). Evolutionary conservation-evaluating the adaptive potential of species. *Evol. Appl.* 7, 963–967. doi: 10.1111/eva.12227

Elshire, R. J., Glaubitz, J. C., Sun, Q., Poland, J. A., Kawamoto, K., Buckler, E. S., et al. (2011). A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. *PLoS ONE* 6:e19379. doi: 10.1371/journal.pone.0019379

Eoche-Bosy, D., Gautier, M., Esquibet, M., Legeai, F., Bretaudreau, A., Bouchez, O., et al. (2017). Genome scans on experimentally evolved populations reveal candidate regions for adaptation to plant resistance in the potato cyst nematode *Globodera pallida*. *Mol. Ecol.* 26, 4700–4711. doi: 10.1111/mec.14240

Eves-van den Akker, S., Lilley, C. J., Reid, A., Pickup, J., Anderson, E., Cock, P. J., et al. (2015). A metagenetic approach to determine the diversity and distribution of cyst nematodes at the level of the country, the field and the individual. *Mol. Ecol.* 24, 5842–5851. doi: 10.1111/mec.13434

Faghihi, J., Donald, P. A., Noel, G., Welacky, T. W., and Ferris, V. R. (2010). Soybean resistance to field populations of *Heterodera glycines* in selected geographic areas. *Plant Health Prog.* doi: 10.1094/PHP-2010-0426-01-RS

Fick, S. E., and Hijmans, R. J. (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. *Int. J. Climatol.* 37, 4302–4315. doi: 10.1002/joc.5086

Foll, M. (2012). *BayeScan v2.1 User Manual*. Available online at: http://cmpg.unibe.ch/software/BayeScan/files/BayeScan2.1_manual.pdf

Foll, M., and Gaggiotti, O. E. (2008). A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. *Genetics* 180, 977–993. doi: 10.1534/genetics.108.092221

Fu, Y.-B., Peterson, G. W., and Dong, Y. (2016). Increasing genome sampling and improving SNP genotyping for genotyping-by-sequencing with new combinations of restriction enzymes. *G3* 6, 845–856. doi: 10.1534/g3.115.025775

Garrison, E., and Marth, G. (2012). Haplotype-Based Variant Detection from Short-Read Sequencing. *arXiv [preprint] arXiv:1207.3907*.

Gautier, M. (2015). Genome-wide scan for adaptive divergence and association with population-specific covariates. *Genetics* 201, 1555–1579. doi: 10.1534/genetics.115.181453

Gautier, M., Gharbi, K., Cezard, T., Foucaud, J., Kerdelhué, C., Pudlo, P., et al. (2013). The effect of RAD allele dropout on the estimation of genetic variation within and between populations. *Mol. Ecol.* 22, 3165–3178. doi: 10.1111/mec.12089

Gendron St-Marseille, A.-F. (2013). Le nématode à Kyste du Soja (Heterodera glycines) : Enjeux des Changements Climatiques sur sa Distribution, sa Reproduction et sur les Probabilités de Synchronisme Avec le Soja (Glycine max) au Québec. Master's Thesis, Université de Sherbrooke, Sherbrooke, QC. Available online at: http://hdl.handle.net/11143/7214

Green, C. D., Greet, D. N., and Jones, F. G. W. (1970). The influence of multiple mating on the reproduction and genetics of *Heterodera rostochiensis* and *H. schachtii*. *Nematologica* 16, 309–326. doi: 10.1163/187529270X00333

Grenier, E., Bossis, M., Fouville, D., Renault, L., and Mugniéry, D. (2001). Molecular approaches to the taxonomic position of Peruvian potato cyst nematodes and gene pool similarities in indigenous and imported populations of *Globodera*. *Heredity* 86, 277–290. doi: 10.1046/j.1365-2540.2001.00826.x

Günther, T., and Coop, G. (2013). Robust identification of local adaptation from allele frequencies. *Genetics* 195, 205–220. doi: 10.1534/genetics.113.152462

Hartl, D. L., and Clark, A. G. (2007). *Principles of Population Genetics*. 4th Edn. Sunderland, Massachusetts Sinauer Associates.

Hassan, M. A., Pham, T. H., Shi, H., and Zheng, J. (2013). Nematodes threats to global food security. *Acta. Agric. Scand. B. Soil Plant Sci.* 63, 420–425. doi: 10.1080/09064710.2013.794858

Hazel, J. R., and Williams, E. E. (1990). The role of alterations in membrane lipid composition in enabling physiological adaptation of organisms to their physical environment. *Prog. Lipid Res.* 29, 167–227. doi: 10.1016/0163-7827(90)90002-3

Hegge, A. H. (1957). Soybean cyst nematode, *Heterodera glycines*, in Missouri. *Plant Dis. Rep.* 41:201.

Hertel, P., Daniel, J., Stegehake, D., Vaupel, H., Kailayangiri, S., Gruel, C., et al. (2013). The ubiquitin-fold modifier 1 (Ufm1) cascade of *Caenorhabditis elegans*. *J. Biol. Chem.* 288, 10661–10671. doi: 10.1074/jbc.M113.458000

Hulme, P. E. (2009). Trade, transport and trouble: managing invasive species pathways in an era of globalization. *J. Appl. Ecol.* 46, 10–18. doi: 10.1111/j.1365-2664.2008.01600.x

Hutchison, D. W., and Templeton, A. R. (1999). Correlation of pairwise genetic and geographic distance measures: inferring the relative influences of gene flow and drift on the distribution of genetic variability. *Evolution* 53, 1898–1914. doi: 10.1111/j.1558-5646.1999.tb04571.x

Kamvar, Z. N., Tabima, J. F., and Grünwald, N. J. (2014). Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. *PeerJ* 2:e281. doi: 10.7717/peerj.281

Kass, R. E., and Raftery, A. E. (1995). Bayes factors. *J. Am. Stat. Assoc.* 90, 773–795. doi: 10.1080/01621459.1995.10476572

Kimura, M., and Weiss, G. H. (1964). The stepping stone model of population structure and the decrease of genetic correlation with distance. *Genetics* 49, 561–576.

Kofler, R., Pandey, R. V., and Schlötterer, C. (2011). PoPoolation2: identifying differentiation between populations using sequencing of pooled DNA samples (pool-Seq). *Bioinformatics* 27, 3435–3436. doi: 10.1093/bioinformatics/btr589

Kondrashov, F. A. (2012). Gene duplication as a mechanism of genomic adaptation to a changing environment. *Proc. R. Soc. B.* 279, 5048–5057. doi: 10.1098/rspb.2012.1108

Kristjansson, G. (2010). 09–59.Plant Health Risk Assessment Unit of the Science Advice Division: Canadian Food Inspection Agency. 45p.

Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler Transform. *Bioinformatics* 25, 1754–1760. doi: 10.1093/bioinformatics/btp324

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. *Bioinformatics* 25, 2078–2079. doi: 10.1093/bioinformatics/btp352

Liu, X. H., Li, J. Q., and Zhang, D. S. (1997). History and status of soybean cyst nematode in China. *Int. J. Nematol.* 7, 18–25.

Lu, F., Lipka, A. E., Glaubitz, J., Elshire, R., Cherney, J. H., Casler, M. D., et al. (2013). Switchgrass genomic diversity, ploidy, and evolution: novel insights from a network-based SNP discovery protocol. *PLoS Genet.* 9:e1003215. doi: 10.1371/journal.pgen.1003215

Lynch, M., Bost, D., Wilson, S., Maruki, T., and Harrison, S. (2014). Population-genetic inference from pooled-sequencing data. *Genome Biol. Evol.* 6, 1210–1218. doi: 10.1093/gbe/evu085

Mathew, F., Markell, S., Jantzi, D., Yan, G., Nelson, B., and Helms, T. (2015). *Soybean cyst nematode. Plant Disease Management*. Fargo, ND: North Dakota State University Extension Service), p. 4.
Mimee, B., Duceppe, M.-O., Véronneau, P.-Y., Lafond-Lapalme, J., Jean, M.,
Belzile, F., et al. (2015). A new method for studying population genetics of cyst
nematodes based on Pool-Seq and genomewide allele frequency analysis. *Mol.
Ecol. Resour.* 15, 1356-1365. doi: 10.1111/1755-0998.12412

Mimee, B., Gagnon, A.-È., Colton-Gagnon, K., and Tremblay, É. (2016). Portrait de la situation du nématode à kyste du soja (*Heterodera glycines*) au Québec (2013–2015). *Phytoprotection* 96, 33–42. doi: 10.7202/1038941ar

Mimee, B., Peng, H., Popovic, V., Yu, Q., Duceppe, M.-O., Tétreault, M.-P., et al. (2014). First report of soybean cyst nematode (*Heterodera glycines* Ichinohe) on soybean in the province of Quebec, Canada. *Plant Dis.* 98:429. doi: 10.1094/PDIS-07-13-0782-PDN

Mitchum, M. G., Wrather, J. A., Heinz, R. D., Shannon, J. G., and
Danekas, G. (2007). Variability in distribution and virulence phenotypes
of *Heterodera glycines* in Missouri during 2005. *Plant Dis.* 91, 1473-1476.
doi: 10.1094/PDIS-91-11-1473

Navon, O., Sul, J. H., Han, B., Conde, L., Bracci, P. M., Riby, J., et al. (2013). Rare variant association testing under low-coverage sequencing. *Genetics* 194, 769-779. doi: 10.1534/genetics.113.150169

Nei, M., Niimura, Y., and Nozawa, M. (2008). The evolution of animal chemosensory receptor gene repertoires: roles of chance and necessity. *Nat. Rev. Genet.* 9, 951-963. doi: 10.1038/nrg2480

Nelson, B. D., and Bradley, C. A. (2003). *Soybean cyst nematode (SCN)*. Fargo, ND: North Dakota State University. Available online at: https://www.ndsu.edu/pubweb/$sim$bernelso/soydiseases/cyst.shtml

Niblack, T. L., Colgrove, A. L., Colgrove, K., and Bond, J. P. (2008). Shift in virulence of soybean cyst nematode is associated with use of resistance from PI 88788. *Plant Health Prog.* doi: 10.1094/PHP-2008-0118-01-RS

Nikolic, A., Kojic, S., Knezevic, S., Krivokapic, Z., Ristanovic, M., and Radojkovic, D. (2011). Structural and functional analysis of SMAD4 gene promoter in malignant pancreatic and colorectal tissues: detection of two novel polymorphic nucleotide repeats. *Cancer Epidemiol.* 35, 265-271. doi: 10.1016/j.canep.2010.10.002

Novak, S. J. (2007). The role of evolution in the invasion process. *PNAS* 104, 3671-3672. doi: 10.1073/pnas.0700224104

Pimentel, D., Zuniga, R., and Morrison, D. (2005). Update on the environmental and economic costs associated with alien-invasive species in the United States. *Ecol. Econ.* 52, 273-288. doi: 10.1016/j.ecolecon.2004.10.002

Plantard, O., Picard, D., Valette, S., Scurrah, M., Grenier, E., and Mugniéry, D. (2008). Origin and genetic diversity of Western European populations of the potato cyst nematode (*Globodera pallida*) inferred from mitochondrial sequences and microsatellite loci. *Mol. Ecol.* 17, 2208-2218. doi: 10.1111/j.1365-294X.2008.03718.x

Pluess, A. R., Frank, A., Heiri, C., Lalagüe, H., Vendramin, G. G., and Oddou-Muratorio, S. (2016). Genome-environment association study suggests local adaptation to climate at the regional scale in *Fagus sylvatica*. *New Phytol.* 210, 589-601. doi: 10.1111/nph.13809

Plummer, M., Best, N., Cowles, K., and Vines, K. (2006). CODA: convergence diagnosis and output analysis for MCMC. *R News* 6, 7–11. Available online at: http://oro.open.ac.uk/22547/

Poland, J. A., Brown, P. J., Sorrells, M. E., and Jannink, J. L. (2012). Development of high-density genetic maps for barley and wheat using a novel two-enzyme genotyping-by-sequencing approach. *PLoS ONE* 7:e32253. doi: 10.1371/journal.pone.0032253

R Core Team (2017). *R: A Language and Environment For Statistical Computing*. Vienna, R Foundation for Statistical Computing. Available online at: http://www.R-project.org/

Revell, L. J. (2017). *Phylogenetic Tools for Comparative Biology (and other things)*. Available online at: https://cran.r-project.org/web/packages/phytools/phytools.pdf

Riggs, R. D. (2004). "History and distribution," in *Biology and Management of Soybean Cyst Nematode, 2nd Edn*, eds D. P. Schmitt, J. A. Wrather,

R. D. Riggs (Marceline, MO: Schmitt and Associates of Marceline),
9-40.

Rousset, F. (1997). Genetic differentiation and estimation of gene flow from *F*-statistics under isolation by distance. *Genetics* 145, 1219-1228.

Savage-Dunn, C., Maduzia, L. L., Zimmerman, C. M., Roberts, A. F., Cohen, S., Tokarz, R., et al. (2003). Genetic screen for small body size mutants in *C. elegans* reveals many TGFβ pathway components. *Genesis* 35, 239–247. doi: 10.1002/gene.10184

Schlötterer, C., Tobler, R., Kofler, R., and Nolte, V. (2014). Sequencing pools of individuals - Mining genome-wide polymorphism data without big funding. *Nat. Rev. Genetics* 15, 749–763. doi: 10.1038/nrg3803

Shurtleff, W., and Aoyagi, A. (2010). *History of Soybeans and Soyfoods in Canada (1831–2010): Extensively Annotated Bibliography and Sourcebook*. Lafayette, CA: Soyinfo Center.

Slack, D. A., Riggs, R. D., and Hamblen, M. L. (1972). The effect of temperature and moisture on the survival of *Heterodera glycines* in the absence of a host. *J. Nematol.* 4, 263-266.

Statistics Canada (2017). *Table 001-0017-Estimated Areas, Yield, Production, Average Farm Price and Total Farm Value of Principal Field Crops, in Imperial Units, Annual, CANSIM (Database)*. [Accessed: November 20, 2017]. Available online at: http://www5.statcan.gc.ca/cansim/a26?lang=eng&id=10017

Tardivel, A., Sonah, H., Belzile, F., and O'Donoghue, L. S. (2014). Rapid identification of alleles at the soybean maturity gene E3 using genotyping by sequencing and a haplotype-based approach. *Plant Genome* 7, 1-9. doi: 10.3835/plantgenome2013.10.0034

Tylka, G. L., and Marett, C. C. (2017). Known distribution of the soybean cyst nematode, *Heterodera glycines*, in the United States and Canada: 1954 to 2017. *Plant Health Prog.* 18, 167-168. doi: 10.1094/PHP-05-17-0031-BR

Van Tassell, C., Smith, T. P., Matukumalli, L. K., Taylor, J. F., Schnabel, R. D., Taylor, L., et al. (2008). SNP discovery and allele frequency estimation by deep sequencing of reduced representation libraries. *Nat. Methods* 5, 247-252. doi: 10.1038/nmeth.1185

Watanabe, N., Ishihara, T., and Ohshima, Y. (2007). Mutants carrying two *sma* mutations are super small in the nematode *C. elegans*. *Genes Cells* 12, 603-609. doi: 10.1111/j.1365-2443.2007.01077.x

Wood, J. F., and Ferkey, D. M. (2016). "GRK roles in C. elegans," in *G Protein-Coupled Receptor Kinases*, eds V. V. Gurevich, and J. J. G. Tesmer (New York, NY: Springer), 283-299.

Wright, S. (1943). Isolation by distance. *Genetics* 28, 114-138.

Yu, Q. (2011). "Soybean cyst nematode (*Heterodera glycines* ichinohe)," in *Soybean Physiology and Biochemistry*, ed H. A. El-Shemy (Rijeka: InTech Europe), 461-474. Available online at: http://cdn.intechopen.com/pdfs-wm/22782.pdf

Zheng, J., Li, Y., and Chen, S. (2006). Characterization of the virulence phenotypes of *Heterodera glycines* in Minnesota. *J. Nematol.* 38, 383-390.

Zheng, J. Z., and Chen, S. Y. (2011). Estimation of virulence type and level of soybean cyst nematode field populations in response to resistant cultivars. *J. Entomol. Nematol.* 3, 37–43. Available online at: http://www.academicjournals.org/journal/JEN/article-abstract/791181810204

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Her Majesty the Queen in Right of Canada, as represented by Agriculture and Agri-Food Canada (AAFC). This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
<figure><img src="image_10.png" /></figure>

Check for
updates

# Real-Time PCR for Diagnosing and Quantifying Co-infection by Two Globally Distributed Fungal Pathogens of Wheat

Araz S. Abdullah¹*, Chala Turo¹, Caroline S. Moffat¹, Francisco J. Lopez-Ruiz¹, Mark R. Gibberd¹, John Hamblin² and Ayalsew Zerihun¹

¹Centre for Crop and Disease Management, School of Molecular and Life Sciences, Curtin University, Bentley, WA, Australia, ²Institute of Agriculture, University of Western Australia, Crawley, WA, Australia

Co-infections – invasions of a host-plant by multiple pathogen species or strains – are common, and are thought to have consequences for pathogen ecology and evolution. Despite their apparent significance, co-infections have received limited attention; in part due to lack of suitable quantitative tools for monitoring of co-infecting pathogens. Here, we report on a duplex real-time PCR assay that simultaneously distinguishes and quantifies co-infections by two globally important fungal pathogens of wheat: *Pyrenophora tritici-repentis* and *Parastagonospora nodorum*. These fungi share common characteristics and host species, creating a challenge for conventional disease diagnosis and subsequent management strategies. The assay uses uniquely assigned fluorogenic probes to quantify fungal biomass as nucleic acid equivalents. The probes provide highly specific target quantification with accurate discrimination against non-target closely related fungal species and host genes. Quantification of the fungal targets is linear over a wide range (5000–0.5 pg DNA µl⁻¹) with high reproducibility (RSD ≤ 10%). In the presence of host DNA in the assay matrix, fungal biomass can be quantified up to a fungal to wheat DNA ratio of 1 to 200. The utility of the method was demonstrated using field samples of a cultivar sensitive to both pathogens. While visual and culture diagnosis suggested the presence of only one of the pathogen species, the assay revealed not only presence of both co-infecting pathogens (hence enabling asymptomatic detection) but also allowed quantification of relative abundances of the pathogens as a function of disease severity. Thus, the assay provides for accurate diagnosis; it is suitable for high-throughput screening of co-infections in epidemiological studies, and for exploring pathogen–pathogen interactions and dynamics, none of which would be possible with conventional approaches.

**Keywords:** co-infections, qPCR, hydrolysis probes, molecular disease diagnosis, pathogen interactions

## OPEN ACCESS

**Edited by:**
Christophe Le May,
Agrocampus Ouest, France

**Reviewed by:**
Martin Chilvers,
Michigan State University,
United States
Syama Chatterton,
Agriculture and Agri-Food Canada
(AAFC), Canada

*Correspondence:
Araz S. Abdullah
araz.abdullah@curtin.edu.au

**Specialty section:**
This article was submitted to Plant Microbe Interactions, a section of the journal Frontiers in Plant Science

**Received:** 23 April 2018
**Accepted:** 05 July 2018
**Published:** 09 August 2018

**Citation:**
Abdullah AS, Turo C, Moffat CS, Lopez-Ruiz FJ, Gibberd MR, Hamblin J and Zerihun A (2018) Real-Time PCR for Diagnosing and Quantifying Co-infection by Two Globally Distributed Fungal Pathogens of Wheat. Front. Plant Sci. 9:1086. doi: 10.3389/fpls.2018.01086

## INTRODUCTION

Co-infections, whereby a host-plant is invaded by multiple pathogens or multiple strains of the same pathogen, are common in the field and can have major consequences for disease ecology and pathogen evolution (Alizon et al., 2013). Despite the recognition of the significance of co-infection (Savary and Zadoks, 1992; Zadoks, 1999; Tack et al., 2012), empirical studies are still few, due mainly to the complexity of distinguishing/quantifying multiple pathogens, which often requires
suitable molecular tools (Tollenaere et al., 2012; Tollenaere et al., 2016). As a result, the theoretical understanding of co-infection has largely outpaced experimental studies in natural and agricultural systems (Hood, 2003; Abdullah et al., 2017) and our capacity to understand and manage the processes and consequences of co-infection is limited until current and future theoretical models can be validated by reliable data.

Methods used for distinguishing multiple pathogens from co-infected tissues involve detection of nucleic acid targets using metagenome sequencing and/or analysis of melt-curves from PCR primers with various annealing properties (Brandfass and Karlovsky, 2006; Tollenaere et al., 2012; Gelaye et al., 2017). Such methods enable multiple pathogen detections but provide limited quantitative information on the relative abundance of each pathogen and their effects on overall disease (Bernreiter, 2017). The quantification of pathogen abundance is critical, as many pathogens are naturally present within plants, but their infection levels, pathogenicity and hence relative impacts can differ vastly (Bakker et al., 2014). At present, real-time quantitative PCR (qPCR) is the most reliable technique for measuring disease load in a sample while providing species specificity (Holland et al., 1991; Bates et al., 2001; Schena et al., 2006). Successful qPCR-based multiplexed assays have been developed for detection of *Phytophthora* diseases of soybean (Rojas et al., 2017), but the accuracy of quantification is often compromised in samples with unbalanced target ratios (Klerks et al., 2004; Atallah et al., 2007). Thus, to ensure reliable quantification, co-infection studies often analyze the abundance of multiple pathogens in separate reactions. Due to the cost of labor and resources required, this approach is limited to small-scale investigations limiting our capacity to investigate co-infections over large-scale host populations.

In order to address the limitations of current approaches, we investigate the use of dual-labeled species-specific probes to simultaneously detect and quantify two globally distributed fungal pathogens of wheat: *Pyrenophora tritici-repentis* and *Parastagonospora nodorum*. These two foliar fungal pathogens cause tan (yellow) spot and septoria nodorum blotch in wheat, respectively, damaging photosynthetically active leaf area and causing substantial yield losses (Bhathal et al., 2003). Symptoms caused by these two diseases are difficult to distinguish and may be misdiagnosed as other unrelated physiological abnormalities (Moffat et al., 2015), creating a challenge for conventional disease diagnosis and subsequent management strategies. Symptoms can also vary markedly depending on host genotype, further hampering accurate disease diagnosis (Lamari and Bernier, 1989). Interestingly, *P. tritici-repentis* carries a pathogenicity/virulence gene, *ToxA*, thought to have been acquired laterally from *Pa. nodorum* (Friesen et al., 2006). This suggests that co-infection by these two fungi is likely to occur in nature and may have consequences for disease management strategies. Despite the high likelihood of co-infection by these two fungi and difficulties of identification, no tools have been developed for the diagnosis of their abundance in co-infected host materials.

Here, we report a duplex assay that distinguishes specific DNA sequences unique to each of these fungi and quantifies their presence by direct comparison to standards amplified in parallel

reactions. The basis of this method is the uniquely assigned fluorogenic reporters for each species sequence. The fluorogenic reporters, when bound to the target, quantify fungal b as nucleic acid equivalents. We report on a series of experiments designed to demonstrate the specificity and sensitivity of the method using both a simulated genomic DNA matrix of both fungi at varying ratios as well as naturally infected wheat leaves collected from the field. We also examine the capacity of the technique to show that a species that is present at very low level can be detected and quantified, with very limited interference, in the presence of an abundant species. Finally, we correlate the amount of fungal DNA quantified using the duplex qPCR assay developed here with the level of infection measured by conventional disease scores.

# MATERIALS AND METHODS

## PCR Primers and Conditions

*P. tritici-repentis* primers were designed to target a species-specific multicopy genomic region described previously (Moffat et al., 2015; See et al., 2016). Briefly, a 4.65-kb region of the *P. tritici-repentis* isolate Pt-1CBFP carrying the *ToxA* gene (supercontig_1.4) was aligned to an orthologous genomic region of the *Pa. nodorum* isolate SN15 (scaffold_55). All known isolates of the targeted fungi were included in this step to ensure that the primers amplify all known isolates of the pathogens of interest. A pair of primers that target a 99-bp fragment located 701-bp upstream of the *ToxA* coding region was designed to detect *P. tritici-repentis* (Table 1). These primers amplify a short fragment (99-bp) within the promoter region of a low molecular weight host-selective toxin. This multicopy region has been detected in a number races of *P. tritici-repentis* isolates collected from around the world (Moolhuijzen et al., 2018). The size of the primer pair was restricted to 99-bp to ensure comparable amplicon sizes between *P. tritici-repentis* and *Pa. nodorum*.

*Pa. nodorum* DNA was amplified using primers modified from those previously described by Oliver et al. (2008). Additional 3-bp was included at the beginning of each primer sequence to allow the probe to overcome the issue of primer-dimer association (Barbisin et al., 2009). *Pa. nodorum* primers amplify a 112-bp fragment of a highly conserved anonymous

<table><caption>TABLE 1 | List of the oligonucleotides used in this study.</caption><thead><tr><th>Sequence ID*</th><th>5'→3' sequence</th><th>Product length (bp)</th><th>GC (%)</th></tr></thead><tbody><tr><td><i>Ptr</i>-Forward</td><td>GTCTCCTCTGGTGGTATG</td><td rowspan="2">99</td><td rowspan="2">55.6</td></tr><tr><td><i>Ptr</i>-Reverse</td><td>GCTCTTAGTGAAGTTCAATC</td></tr><tr><td><i>Ptr</i>-Probe</td><td>TACCTCTACTCGTCGCCTATGG</td><td rowspan="3">112</td><td rowspan="3">45.8</td></tr><tr><td><i>Pn</i>-Forward</td><td><b>ACCGCATTACCAAAGTTC</b></td></tr><tr><td><i>Pn</i>-Reverse</td><td><b>ACTGGAACTGGAACAATAAG</b></td></tr><tr><td><i>Pn</i>-Probe</td><td>CCTGAATGCTCTTGACACTTGGTT</td><td></td><td></td></tr></tbody></table>

<footnote>*Ptr and Pn refer to *P. tritici-repentis* and *Pa. nodorum*, respectively. Sequences in bold are pre-published primers modified from Oliver et al. (2008). Other sequences are designed in this study. Primers and probes were designed to amplify short fragments of the targeted pathogen DNA.</footnote>
gene (SNOG_01116-1). This gene has no significant similarity to any other sequences in the publicly available genome databases (Oliver et al., 2008). All primers were designed using OligoArchitect™ primer analyzer (Sigma, Life Science) and scanned against the National Center for Biotechnology Information GenBank database using basic local alignment search tools to ensure their specificity.

The designed primers were subjected to conventional PCR to confirm their specificity. Each PCR reaction contained 1xMyTaq buffer, 250 nM forward primer, 250 nM reverse primer, 1 unit MyTaq DNA polymerase (Bioline), and 5 ng DNA template. Reactions were performed as follows: 3 min initial denaturation at 95°C, 35 cycles of 30 s denaturation at 95°C, 30 s of annealing at 58°C and 1 min extension at 72°C. Electrophoresis of PCR products was performed on 2% agarose gels stained with SybrSafe (Life Technologies) and visualized under UV light. Reproducibility of the results was confirmed by running the PCR with negative controls in duplicates. The PCR step also included DNA samples from five common fungal pathogens of cereals as negative controls. Colonies of *Pyrenophora teres* f. sp. *maculata*, *Pyrenophora teres* f. sp. *teres*, *Blumeria graminis* f. sp. *tritici*, *Alternaria alternata*, and *Fusarium graminearum* were kindly provided by Steven Chang, Curtin University/Centre for Crop and Disease Management. Species identity was confirmed by sequencing the internal transcribed spacer of the ribosomal DNA (Schoch et al., 2012). All sequence analyses and multiple sequence alignments were carried out using Geneious version R6.1.6.

## Real-Time PCR Probes and Conditions

Two dual-labeled probes were custom-designed and assigned to hybridize with a complementary region between the forward and the reverse primers. The *P. tritici-repentis* and *Pa. nodorum* probes were 23-bp and 20-bp in length, respectively. The length of each probe was chosen to ensure probe-primer hybrids are formed in a complementary manner with the length of the primers. Each probe was labeled with a unique fluorogenic reporter to ensure that target sequences of both pathogens were amplified simultaneously but detected independently. The *P. tritici-repentis* probe was labeled with 6-carboxyfluorescein (FAM™; Sigma-Aldrich). The *Pa. nodorum* probe was labeled with CAL Fluor Gold® (CFG; Sigma-Aldrich). FAM has emission maxima between 494 nm and 518 nm and CFG emission peaks between 538 nm and 559 nm. The fluorogenic reporters were selected based on the capacity of the CFX96 detection system, the instrument used in this study, to resolve overlapping spectra. This was determined prior to carrying out the experiments using a spectra overlay tool available at http://www.Qpcrdesign.com/spectral-overlay. Both probes were paired with the non-fluorescent black hole quencher-1 (BHQ-1®; Sigma-Aldrich).

Probes and their matching primers were run on a 96-well spectrofluorometric thermal cycler (Bio-Rad CFX96) with the following conditions: 15 min at 95°C, 15 s denaturation at 95°C, 20 s at 72°C followed by 40 cycles of 15 s at 95°C and 30 s at 58°C. Each 20 µl reaction volume contained 5 µl of the sample and 12 µl iQ™ Multiplex Powermix (Bio-Rad). For a fixed amount of target template, *P. tritici-repentis* DNA was amplified faster than *Pa. nodorum* DNA. Hence, probe and

primers of *P. tritici-repentis* were restricted to obtain comparable quantification cycles (Cq) to that of *Pa. nodorum*. In a preliminary experiment, the concentrations of primers and probes were gradually increased from 50 up to 450 nM at 50 nM intervals. *P. tritici-repentis* DNA was amplified using 200 nM forward primer, 200 nM reverse primer and 100 nM probe. *Pa. nodorum* DNA was amplified using 250 nM forward primer, 250 nM reverse primer and 150 nM probe. Unless specified, reactions were carried out in duplex where primers and probes for both species were applied together. A preliminary experiment showed a comparable amplification efficiency between duplex and singleplex reactions and no evidence of cross-amplification among primers and probes of the two species was observed (Supplementary Table 1). Presence of any non-specific amplicon was examined using post-PCR melt curve analysis.

## DNA Extraction and Quantification

Pure genomic DNA from fungal colonies was extracted using the Bio-sprint 15 plant DNA kit (Qiagen) as per the manufacturer's protocol. Fungal cultures were maintained on agar plates as described elsewhere (Moffat et al., 2015). Mycelia were harvested from these cultures and ground into a fine powder in liquid nitrogen. Subsamples (40 mg ground tissue) were placed into 1.5-ml microtubes, which were then used for DNA extraction. The concentration of DNA in each subsample was determined using a NanoDrop 2000 UV-Vis spectrophotometer (Thermo Scientific) and diluted to 50 ng µl⁻¹ in ultrapure PCR grade water. DNA was stored at -20°C until used. Aerosol protected pipette tips were used throughout the extraction and quantification steps to prevent DNA contamination. DNA from infected wheat leaves was extracted using the same DNA extraction kit and quantified as described above.

## DNA Spiking and Field Validation

In two independent experiments, a fixed concentration of fungal DNA of one pathogen (5 ng µl⁻¹) was spiked into a progressively decreasing DNA concentration of the other pathogen (5, 0.5, 0.05, and 0.005 ng DNA µl⁻¹). This generated a ratio of DNA concentration from one pathogen to the other in the sample ranging from 1:1 down to 1:10000. The spiking aimed to simulate the analysis of samples derived from infection conditions whereby the two species occur at different relative abundance, which is typical of many foliar fungal infections. Recovery of fungal DNA for the respective pathogens was expressed as the ratio of the total concentration of fungal DNA quantified using the qPCR method to the amount of DNA added × 100.

In a third experiment, recovery of fungal DNA was measured in fungal DNA samples spiked with abundant wheat DNA. The concentration of fungal DNA in the samples was progressively decreased from 5 down to 0.05 ng µl⁻¹ while background wheat DNA was increased from 5 to 100 ng µl⁻¹. The experiment simulated analytical conditions where fungal DNA is present in small concentrations against a background of ample wheat DNA, such as would occur where infection severity and/or incidence are low.

To evaluate how well the assay works for field samples, diseased leaves from the wheat variety Scout were collected from
a site in the southwest of Western Australia (31°.74S, 116°.70E). Scout is rated as susceptible to very susceptible to both *P. tritici-repentis* and *Pa. nodorum* (DAFWA, 2016). Sampled leaves were visually inspected for diseased leaf area and given scores on 0-100 scale. Leaves were then split into two groups; one group (n = 9) was surface-sterilized in 2% chlorine and incubated on agar Petri-dishes in an attempt to characterize the causal agent of the disease. Leaves from the second group (n = 9) were ground in liquid nitrogen and used for DNA extraction. 50 ng µl⁻¹ gDNA from the infected leaves, along with 50 ng µl⁻¹ gDNA from uninfected leaves from glasshouse-grown plants (negative controls), were first analyzed by conventional PCR. A further 50 ng µl⁻¹ gDNA from the same infected leaves and controls were then analyzed using qPCR.

## Detection and Quantification of Fungal Biomass

Fluorescence data from the qPCR machine were retrieved during the annealing step of every Cq. Threshold fluorescence was set automatically by the instrument manager system (CFX Manager Version 2.0.8) before carrying out the assay. A log-linear standard curve of a 10-fold dilution series corresponding to 5000 to 0.5 pg µl⁻¹ was generated by plotting the logarithms of known concentrations of fungal DNA against the Cq values. The resulting regression equations were used to calculate fungal DNA in unknown samples. No-template controls, where sterile water was added instead of DNA, were included in each reaction. Limit of detection of the duplex assay, the lowest concentration at which reliable detection can be achieved, was determined following Armbruster and Pry (2008). All reactions were run with three replicates and samples that gave positive fluorescence before no-template controls were considered positive.

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 1** | Specificity testing of primers by agarose gel electrophoresis. (A) Primers of *P. tritici-repentis*, *Ptr-Forward* and *Ptr-Reverse*, were tested against *Pa. nodorum* DNA, and (B) primers of *Pa. nodorum*, *Pn-Forward* and *Pn-Reverse*, were tested against *P. tritici-repentis* DNA in two separate PCR reactions. Reactions included negative controls of DNA samples from five cereal fungal pathogens. Reactions including positive and negative controls were electrophoresed on 2% agarose with two technical replicates.</figcaption></figure>

## RESULTS

### Specificity of the Assay

To test how well *P. tritici-repentis* and *Pa. nodorum* could be distinguished from each other, as well as from other common wheat pathogens and the host, we carried out two conventional PCR assays each using 5 ng gDNA. Five non-target controls of DNA from closely related cereal fungal pathogens and DNA from the wheat cultivar Scout were included in this step. Reactions were run in separate wells (i.e., singleplex). Conventional PCR provided amplicons of the expected sizes, and *P. tritici-repentis* and *Pa. nodorum* were distinguished based on the size of the amplicons (Supplementary Figure 1). *P. tritici-repentis* and *Pa. nodorum* primers only amplified the respective pathogen DNA. None of the five non-target controls or DNA from wheat gave specific amplicons that could be detected by conventional PCR (Figures 1A,B).

We then tested the specificity of the probes and their matching primers in two-singleplex real-time qPCR assays each containing 5 ng gDNA. Detection only occurred for probes that were complementary to the expected sequences, and none of the five non-target controls or the wheat DNA

had specific amplification during 40 real-time qPCR cycles (Figures 2A,B). Samples that included no-template DNA or those that contained wheat DNA were negative during the course of the reaction (Figures 2A,B). *P. tritici-repentis* and *Pa. nodorum* were distinguished based on the emission spectra of the fluorogenic reporters that were used to label each pathogen probe.

### Dynamic Range, Efficiency, and Reproducibility of the Assay

We carried out two duplexed standard curve experiments on DNA of *P. tritici-repentis* and *Pa. nodorum* mixed at equal ratios. A 10-fold dilution series in the range of 5000 to 0.5 pg µl⁻¹ was generated, analyzed and plotted against the number of Cq required to detect fluorescence signals. Each dilution was prepared with three replicates except the no-template control, which was analyzed using 10 replicates. Fitting log-linear standard curves between Cq and fungal DNA resulted in correlation coefficients (r) -0.997 and -0.996 for *P. tritici-repentis* and *Pa. nodorum*, respectively. The standard curves for both species were linear over 10,000-fold dilution, and the calculated amplification efficiency was 91.28% for *P. tritici-repentis* and 99.25% for *Pa. nodorum* (Figures 3A,B). The slopes and intercepts of the log-linear curves were comparable between *P. tritici-repentis* and *Pa. nodorum* (p > 0.05). Mean Cq value (±standard deviation bars, where visible) of three replicated reactions was highly reproducible with relative standard deviation ≤ 10%
<figure><img src="image_13.png" /><figcaption>**FIGURE 2** | Specificity testing of primers and their matching probes in singleplex real-time qPCR settings. Amplification curve for *P. tritici-repentis* is plotted in blue (*A*) and amplification curve of *Pa. nodorum* is plotted in orange (*B*). Primers and probes (**Table 1**) were tested against DNA from the five fungal species, a negative wheat control, and a no template control and reactions were run separately. Data represent means ± standard deviation (n = 3).</figcaption></figure>

(Figures 3A,B). The variation around the mean was independent of template concentration, and the Cq value for 10 no-template controls averaged 36.618 with a relative standard deviation of 0.35%.

## Sensitivity and Limit of Detection of the Assay

We carried out two spiking experiments using a simulated matrix of *P. tritici-repentis* and *Pa. nodorum* DNA mixed at various ratios. Results of the spiking experiments demonstrated that reducing DNA ratio of one species while keeping DNA

of the other species constant did not affect the detection limit of the assay. We were able to accurately quantify fungal DNA of either pathogen down to a relative ratio of 1:100 (Figures 4A,B). However, further reduction of the DNA ratio to 1:1000 underestimated the amount of *P. tritici-repentis* by 11.75% and *Pa. nodorum* by 18.55%. This underestimation became more pronounced with further dilution, and at 100 ppm (i.e., DNA ratios of 1:10000), neither pathogen was quantifiable in the presence of the other (Figures 4A,B). The minimal DNA concentration for detection was 0.059 pg µl⁻¹ for *P. tritici-repentis* and 0.036 pg µl⁻¹ for *Pa. nodorum*.
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**FIGURE 3** | (A) The relationship between quantification cycle and logarithm of the concentration of fungal DNA in duplexed qPCR settings. Triplicate dilution series corresponding to gDNA concentrations of 5000, 500, 50, 5, and 0.5 pg µl⁻¹ were prepared. No-template samples were included in every reaction as negative controls (n = 10). The quantification cycle at which fluorescent signals were observed is plotted against the logarithm of DNA concentrations of *P. tritici-repentis* (A) and *Pa. nodorum* (B). The corresponding regression equations and coefficient of determinations (R²) are shown on the plot. Data are means ± standard deviation where visible (n = 3).</figcaption></figure>

We also tested the ability of the assay to detect and quantify
fungal DNA in a sample that included abundant wheat DNA.
We were able to quantify both species with high accuracy up
to a *P. tritici-repentis*: *Pa. nodorum*: wheat DNA ratio of 1:1:200
(Figure 5). Recovery of fungal DNA, total amount of DNA
quantified using the qPCR method, was high at 95.69% for
*P. tritici-repentis* and 94.20% for *Pa. nodorum* when the fungal to
wheat DNA ratios dropped to 1:1:2000. At these ratios, however,
relative standard deviations were 10.23% and 13.90% for *P. tritici-
repentis* and *Pa. nodorum*, respectively, affecting the reliability of
fungal DNA quantification.

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**FIGURE 4** | Detection and quantification of P. tritici-repentis and *Pa. nodorum* in a simulated DNA matrix of various ratios. (A) *P. tritici-repentis* DNA was spiked with *Pa. nodorum* DNA at ratios of 1:1, 1:10, 1:100, 1:1000, and 1:10000. (B) *Pa. nodorum* was spiked with *P. tritici-repentis* DNA at same ratios. The starting concentration was 5 ng µl⁻¹ of each pathogen DNA. The subsequent ratios were sequential 10-fold dilutions (0.5, 0.05, 0.005, and 0.0005 ng µl⁻¹) of one pathogen DNA, with the other held constant at 5 ng µl⁻¹. Data represent means ± standard deviation (n = 6).</figcaption></figure>

Field Evaluation of Naturally Infected
Plants

We collected eighteen naturally infected wheat leaves from a
site in the southwest of Western Australia. These leaves showed
typical symptoms of tan spot with distinct yellow halos and tan
chlorotic lesions (Figure 6A). No apparent symptoms of *Pa.
nodorum* were visible on these leaves. Leaves were split into two
groups; one group (n = 9) was incubated on agar Petri-dishes
in an attempt to characterize the causal agent of the disease.
Incubation on agar only yielded colonies typical of *P. tritici-
repentis* (data not shown).

DNA was extracted from the second leaf group (n = 9) and analyzed using conventional PCR. Only one out of nine duplicated PCR reactions tested positive to *Pa. nodorum*. All reactions were positive to *P. tritici-repentis*. DNA from the same leaves that were used in the conventional PCR was then analyzed using qPCR. The duplex qPCR assay was able to detect signals corresponding to sequences of *P. tritici-repentis* and *Pa. nodorum* despite the absence of any visible symptoms of
<figure><img src="image_18.png" /></figure>

**FIGURE 5** | Simulated DNA matrix of various ratios of fungi to wheat DNA. Increasing concentrations (5–100 ng µl⁻¹) of wheat DNA were spiked with decreasing concentrations of fungal DNA (5–0.05 ng µl⁻¹). Quantification was done in real-time quantitative PCR and reactions were duplexed each containing four fungal-specific primers and two complementary probes (Table 1). DNA concentrations respectively were 5:5:5, 5:5:10, 5:5:100, 0.5:0.5:100, and 0.05:0.05:100 ng µl⁻¹. Data represent means ± standard deviation (n = 6). Wheat DNA was not quantified/detected.

*Pa. nodorum* on these leaves (Figure 6B). There was a good agreement (r = -0.82) between total fungal DNA quantified by qPCR and conventional disease scores that were collected in the field (Figure 6B; black circles/line). Conventional disease scores were significantly correlated with the increase in *P. tritici-repentis* DNA, and *P. tritici-repentis* DNA was greater contributor to the disease scores than *Pa. nodorum* (Figure 6B; blue and orange circles/lines).

DISCUSSION

Quantification of the relative abundance of co-infecting pathogens requires selection of primers and probes that are compatible with each other whilst distinguishing between the co-infecting species. Recent work has identified a 235-bp multicopy region present in the *P. tritici-repentis* genome, and primers designed in this region were able to detect *P. tritici-repentis* from wheat leaves, even prior to the visible appearance of tan spot lesions (See et al., 2016). However, these primers cross-hybridized with *Pyrenophora teres* f. sp. teres, which the authors suggested may be overcome by the addition of fluorescence-labeled probes in the middle of the amplicon.

In this study, primers were chosen to produce relatively smaller amplicons (99–112 bp) that are more efficient for amplifying specific products than longer primers (Bustin et al., 2009). Our primers successfully distinguished targeted pathogens from non-target closely related fungal species and host genes. However, for increased specificity, fluorogenic probes were

designed to hybridize to sequences within the primers enabling the assay to simultaneously distinguish and quantify DNA associated with *P. tritici-repentis* from that associated with *Pa. nodorum*. The length of the probes was designed to ensure that stable primer-probe hybrids were formed at the same annealing temperature for both species, irrespective of the length of the primers. Furthermore, probes were labeled with a unique fluorogenic-reporter allowing DNA sequences from *P. tritici-repentis* and *Pa. nodorum* to be amplified simultaneously but independent of each other. The choice of the fluorogenic reporter was decided based on instrument capability in resolving overlapping spectra. FAM-labeled probes are excited at lower emission spectra than CFG-labeled probes and hence are expected to produce stronger signals. Nevertheless, restricting the concentration of the FAM-labeled probe in the reaction optimized the fraction of the amplicon that is bound to the probe and produced comparable signals for *P. tritici-repentis* and *Pa. nodorum*. In addition to the fluorogenic reporters, probes were labeled with dark-hole quenchers to inhibit probe signals when probes are free in the solution (Holland et al., 1991). The dark-hole quencher was sufficient to inhibit signals from probes that are not in perfect contact with the target, resulting in accurate species-specific quantifications.

A series of experiments were conducted demonstrating the specificity, efficiency, and reproducibility of this assay for simultaneous detection and quantification of *P. tritici-repentis* and *Pa. nodorum*. The method reported here resolved two signals each unique to the assayed fungi. No cross-amplification with either host plant DNA or DNA from other closely related fungal
A

<figure><img src="image_19.png" /></figure>

B

<figure><img src="image_20.png" /><figcaption>**FIGURE 6** | Quantification of fungal DNA in naturally infected wheat leaves. (A) Leaves displaying tan spot symptoms with tan necrotic centers and yellow halos. (B) A linear model fitted into the relationship between fungal DNA measured in real-time quantitative PCR and conventional disease score. The corresponding regression equation and coefficient of determination (R²) are shown on the plot. ns and ** refer to not significant and significant (*P* < 0.01), respectively.</figcaption></figure>

pathogens of cereals was observed. The assay amplified the targeted sequences with high specificity, and we were able to detect the presence of *P. tritici-repentis* at 0.059 pg µl⁻¹ and *Pa. nodorum* at 0.036 pg µl⁻¹ in mixed fungal DNA matrices. Furthermore, in mixed gDNA matrices of both targeted species, the assay was able to detect fungal signals at a relative abundance as low as 1 in 10,000 with ≤10% RSD. Similarly, the assay was able to quantify the presence of *Pa. nodorum* DNA down to 23-pg in leaf samples naturally co-infected by abundant *P. tritici-repentis* up to 2192-pg. Nonetheless, there was a slight interference to the quantification of the less abundant target when DNA of the other target and/or wheat DNA was highly abundant. This may occur because reaction components including DNA polymerase, dNTPs, and MgCl₂ become limiting in later qPCR cycles. The more abundant target may compete for reaction components with the less abundant target, delaying its amplification (Bustin et al., 2009). As a result, the quantification of the less abundant target, although remaining specific

and within the detection limit, may be compromised. One way to minimize this interference is to optimize concentrations of reaction components sequentially (Bustin et al., 2009). Regardless, the assay is highly stable and reproducible. The reproducibility was confirmed by the small standard deviations of triplicates samples in the standard curves (≤10%). The amplification curves of replicate runs showed strong overlap indicating stability and reproducibility.

A limitation of most nucleic acid-based quantification assays is an inability to distinguish between DNA from living and dead cells (Pinheiro et al., 2016). Although fungal cell viability was not assessed in this work, the high correlation coefficient, 0.826, especially between *P. tritici-repentis* DNA and the conventional disease score, suggests that viability may not have a significant impact on measuring disease as DNA equivalent. Quantification of fungal DNA has been shown to provide accurate measurement of disease severity for *P. tritici-repentis* (See et al., 2016) and *Pa. nodorum* (Oliver et al., 2008). Another limitation more specific to qPCR is that only a few targets can readily be accommodated in a single homogeneous assay. The number of targets that can be fitted in one assay depends on the design of the assay and the ability of the instrument to resolve overlapping spectra. The instrument used in this study has a blue monochromatic laser source for the excitation of fluorophores, and can simultaneously detect five targets on five different channels. This limits the fluorophores that could be efficiently excited by the blue laser to those that are distant from the far-red wavelengths. Instruments with high resolution-detectors and variable-wavelength light sources may allow more targets to be included in an assay. Instruments with such specifications are available and can be employed for detection and quantification of large-scale screening of co-infection by multiple pathogens.

Detection of multiple pathogens has also been achieved by the use of post-PCR dissociation curve analysis. Such methods use a set of primers with different melting temperatures and G + C contents. Using this technique up to 10 microbial pathogens of humans have been distinguished in a single assay (Gelaye et al., 2017; Jiang et al., 2017). Post-PCR methods, although specific and accommodating several targets in a single assay, do not provide quantitative information on the level of disease. Plant tissues can be inhabited by a large number of pathogenic and non-pathogenic microbes. Nevertheless, infection levels of these microbes and the damage they cause to the plant can vary greatly (Bakker et al., 2014). Our assay offers several advantages over the dissociation curve analysis and many other methods for diagnosis of co-infection. First, no post-reaction processing is required. Reactions are driven to completion allowing higher levels of sensitivity even when target sequences are present at low concentrations. Second, the assay provides quantitative data on the level of disease load in a sample, which other methods including metagenome sequencing do not offer. Finally, the cost-effectiveness, although not directly estimated, of testing multiple agents in a single test allows more testing using the same amount of reagents and staff time. When coupled with a 96-well capacity, this assay offers a sensitive and quantitative high throughput methodology for detection of co-infection in plants. The method
can also be used to study epidemiological consequences of co-infection in the field.

## AUTHOR CONTRIBUTIONS

AA contributed to the conception of the idea, defined/selected the primers, conducted the experiments, analyzed the data, and drafted the manuscript. CT assisted in setting up the conventional PCR and testing of the primers. CM, FL-R, AZ, MG, and JH supervised the work and provided critical suggestions on the article. All authors read and approved the final version of the article.

## FUNDING

AA is a recipient of Research Training Program (RTP), an Australian government funded scholarship for Ph.D. students. The Centre for Crop and Disease Management (CCDM) defrayed the publishing costs for this review. CCDM is a national research

## REFERENCES

* Abdullah, A. S. A., Moffat, C., Lopez Ruiz, F. J. L., Gibberd, M. G., Hamblin, J., and Zerihun, A. Z. (2017). Host-multi-pathogen warfare: pathogen interactions in co-infected plants. *Front. Plant Sci.* 8:1806. doi: 10.3389/fpls.2017.01806

* Alizon, S., De Roode, J. C., and Michalakis, Y. (2013). Multiple infections and the evolution of virulence. *Ecol. Lett.* 16, 556-567. doi: 10.1111/ele.12076

* Armbruster, D. A., and Pry, T. (2008). Limit of blank, limit of detection and limit of quantitation. *Clin. Biochem. Rev.* 29(Suppl. 1), S49-S52.

* Atallah, Z. K., Bae, J., Jansky, S. H., Rouse, D. I., and Stevenson, W. R. (2007). Multiplex real-time quantitative PCR to detect and quantify *Verticillium dahliae* colonization in potato lines that differ in response to *Verticillium wilt*. *Phytopathology* 97, 865-872. doi: 10.1094/PHYTO-97-7-0865

* Bakker, M. G., Schlatter, D. C., Otto-Hanson, L., and Kinkel, L. L. (2014). Diffuse symbioses: roles of plant-plant, plant-microbe and microbe-microbe interactions in structuring the soil microbiome. *Mol. Ecol.* 23, 1571-1583. doi: 10.1111/mec.12571

* Barbisin, M., Fang, R., O'shea, C. E., Calandro, L. M., Furtado, M. R., and Shewale, J. G. (2009). Developmental validation of the quantifiler® duo DNA quantification kit for simultaneous quantification of total human and human male DNA and detection of PCR inhibitors in biological samples. *J. Forensic Sci.* 54, 305-319. doi: 10.1111/j.1556-4029.2008.00951.x

* Bates, J., Taylor, E., Kenyon, D., and Thomas, J. (2001). The application of real-time PCR to the identification, detection and quantification of *Pyrenophora* species in barley seed. *Mol. Plant Pathol.* 2, 49-57. doi: 10.1046/j.1364-3703.2001.00049.x

* Bernreiter, A. (2017). Molecular diagnostics to identify fungal plant pathogens-a review of current methods. *Ecuad. Calid. Rev. Cient. Ecuat.* 4, 26-35.

* Bhathal, J., Loughman, R., and Speijers, J. (2003). Yield reduction in wheat in relation to leaf disease from yellow (tan) spot and *Septoria nodorum* blotch. *Eur. J. Plant Pathol.* 109, 435-443. doi: 10.1023/A:1024277420773

* Brandfass, C., and Karlovsky, P. (2006). Simultaneous detection of *Fusarium culmorum* and *F. graminearum* in plant material by duplex PCR with melting curve analysis. *BMC Microbiol.* 6:4. doi: 10.1186/1471-2180-6-4

* Bustin, S. A., Benes, V., Garson, J. A., Hellemans, J., Huggett, J., Kubista, M., et al. (2009). The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. *Clin. Chem.* 55, 611-622. doi: 10.1373/clinchem.2008.112797

* DAFWA (2016). Wheat Variety Guide for Western Australia. Department of Agriculture and Food Western Australia. Available at: www.agric.wa.gov.au/cropdisease

centre co-funded by Curtin University and the Grains Research and Development Corporation (project grant CUR00023). The funders had no role in ideas, design, decision to publish, or preparation of the manuscript.

## ACKNOWLEDGMENTS

We thank Richard Oliver, Centre for Crop and Disease Management, for critically reading the manuscript and contributing to the genesis of the research. Belinda Cox and Weiwei Deng provided methodological and technical assistance, and Steven Chang shared some of the fungal isolates used in this study.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.01086/full#supplementary-material

* Friesen, T. L., Stukenbrock, E. H., Liu, Z., Meinhardt, S., Ling, H., Faris, J. D., et al. (2006). Emergence of a new disease as a result of interspecific virulence gene transfer. *Nat. Genet.* 38, 953-956. doi: 10.1038/ng1839

* Gelaye, E., Mach, L., Kolodziejek, J., Grabherr, R., Loitsch, A., Achenbach, J. E., et al. (2017). A novel HRM assay for the simultaneous detection and differentiation of eight poxviruses of medical and veterinary importance. *Sci. Rep.* 7:42892. doi: 10.1038/srep42892

* Holland, P. M., Abramson, R. D., Watson, R., and Gelfand, D. H. (1991). Detection of specific polymerase chain reaction product by utilizing the 5'-3' exonuclease activity of *Thermus aquaticus* DNA polymerase. *Proc. Natl. Acad. Sci. U.S.A.* 88, 7276-7280. doi: 10.1073/pnas.88.16.7276

* Hood, M. (2003). Dynamics of multiple infection and within-host competition by the anther-smut pathogen. *Am. Nat.* 162, 122-133. doi: 10.1086/375539

* Jiang, Y., He, L., Wu, P., Shi, X., Jiang, M., Li, Y., et al. (2017). Simultaneous identification of ten bacterial pathogens using the multiplex ligation reaction based on the probe melting curve analysis. *Sci. Rep.* 7:5902. doi: 10.1038/s41598-017-06348-z

* Klerks, M., Zijlstra, C., and Van Bruggen, A. (2004). Comparison of real-time PCR methods for detection of *Salmonella enterica* and *Escherichia coli* O157: H7, and introduction of a general internal amplification control. *J. Microbiol. Methods* 59, 337-349. doi: 10.1016/j.mimet.2004.07.011

* Lamari, L., and Bernier, C. (1989). Evaluation of wheat lines and cultivars to tan spot [*Pyrenophora tritici-repentis*] based on lesion type. *Can. J. Plant Pathol.* 11, 49-56. doi: 10.1080/07060668909501146

* Moffat, C., See, P. T., and Oliver, R. (2015). Leaf yellowing of the wheat cultivar mace in the absence of yellow spot disease. *Australas. Plant Pathol.* 44, 161-166. doi: 10.1007/s13313-014-0335-2

* Moolhuijzen, P., See, P. T., Hane, J. K., Shi, G., Liu, Z., Oliver, R. P., et al. (2018). Comparative genomics of the wheat fungal pathogen *Pyrenophora tritici-repentis* reveals chromosomal variations and genome plasticity. *BMC Genomics* 19:279. doi: 10.1186/s12864-018-4680-3

* Oliver, R. P., Rybak, K., Shankar, M., Loughman, R., Harry, N., and Solomon, P. (2008). Quantitative disease resistance assessment by real-time PCR using the *Stagonospora nodorum*-wheat pathosystem as a model. *Plant Pathol.* 57, 527-532. doi: 10.1111/j.1365-3059.2007.01787.x

* Pinheiro, E. T., Neves, V. D., Reis, C. C., Longo, P. L., and Mayer, M. P. (2016). Evaluation of the propidium monoazide-quantitative polymerase chain reaction method for the detection of viable *Enterococcus faecalis*. *J. Endod.* 42, 1089-1092. doi: 10.1016/j.joen.2016.04.003

* Rojas, J. A., Miles, T. D., Coffey, M. D., Martin, F. N., and Chilvers, M. I. (2017). Development and application of qPCR and RPA genus-and
species-specific detection of *Phytophthora sojae* and *P. sansomeana* root rot pathogens of soybean. *Plant Dis.* 101, 1171-1181. doi: 10.1094/PDIS-09-16-1225-RE

Savary, S., and Zadoks, J. (1992). Analysis of crop loss in the multiple pathosystem groundnut-rust-late leaf spot. III. Correspondence analyses. *Crop Prot.* 11, 229-239. doi: 10.1016/0261-2194(92)90042-4

Schena, L., Hughes, K. J., and Cooke, D. E. (2006). Detection and quantification of *Phytophthora ramorum*, *P. kernoviae*, *P. citricola* and *P. quercina* in symptomatic leaves by multiplex real-time PCR. *Mol. Plant Pathol.* 7, 365-379. doi: 10.1111/j.1364-3703.2006.00345.x

Schoch, C. L., Seifert, K. A., Huhndorf, S., Robert, V., Spouge, J. L., Levesque, C. A., et al. (2012). Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for Fungi. *Proc. Natl. Acad. Sci. U.S.A.* 109, 6241-6246. doi: 10.1073/pnas.1117018109

See, P. T., Moffat, C. S., Morina, J., and Oliver, R. P. (2016). Evaluation of a multilocus indel DNA region for the detection of the wheat tan spot pathogen *Pyrenophora tritici-repentis*. *Plant Dis.* 100, 2215-2225. doi: 10.1094/PDIS-03-16-0262-RE

Tack, A. J., Thrall, P. H., Barrett, L. G., Burdon, J. J., and Laine, A. L. (2012). Variation in infectivity and aggressiveness in space and time in wild host-pathogen systems: causes and consequences. *J. Evol. Biol.* 25, 1918-1936. doi: 10.1111/j.1420-9101.2012.02588.x

Tollenaere, C., Susi, H., and Laine, A.-L. (2016). Evolutionary and epidemiological implications of multiple infection in plants. *Trends Plant Sci.* 21, 80-90. doi: 10.1016/j.tplants.2015.10.014

Tollenaere, C., Susi, H., Nokso-Koivisto, J., Koskinen, P., Tack, A., Auvinen, P., et al. (2012). SNP design from 454 sequencing of *Podosphaera plantaginis* transcriptome reveals a genetically diverse pathogen metapopulation with high levels of mixed-genotype infection. *PLoS One* 7:e52492. doi: 10.1371/journal.pone.0052492

Zadoks, J. (1999). Reflections on space, time, and diversity. *Annu. Rev. Phytopathol.* 37, 1-17. doi: 10.1146/annurev.phyto.37.1.1

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Abdullah, Turo, Moffat, Lopez-Ruiz, Gibberd, Hamblin and Zerihun. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
<figure><img src="image_21.png" /></figure>

Check for
updates

# Defense Compounds Rather Than Nutrient Availability Shape Aggressiveness Trait Variation Along a Leaf Maturity Gradient in a Biotrophic Plant Pathogen

Agathe Maupetit¹, Romain Larbat², Michaël Pernaci¹†, Axelle Andrieux¹, Cécile Guinet³, Anne-Laure Boutigny³†, Bénédicte Fabre¹, Pascal Frey¹ and Fabien Halkett¹*

¹ Université de Lorraine, INRA, IAM, Nancy, France, ² Université de Lorraine, INRA, LAE, Nancy, France, ³ ANSES, Laboratoire de la Santé des Végétaux, Unité de Mycologie, Malzéville, France

**OPEN ACCESS**

**Edited by:**

Christophe Le May,
Agrocampus Ouest, France

**Reviewed by:**

Sabine Banniza,
University of Saskatchewan, Canada
Fournet Sylvain,
INRA (Paris), France

**\*Correspondence:**

Fabien Halkett
fabien.halkett@inra.fr

**†Present address:**

Michaël Pernaci,
PalmElit SAS, Montferriez-sur-Lez,
France
Anne-Laure Boutigny,
ANSES – Laboratoire de la Santé des
Végétaux, Unité de Bactériologie,
Virologie et OGM, Angers, France

**Specialty section:**

This article was submitted to
Plant Microbe Interactions,
a section of the journal
Frontiers in Plant Science

Received: 20 April 2018

Accepted: 03 September 2018

Published: 28 September 2018

**Citation:**

Maupetit A, Larbat R, Pernaci M, Andrieux A, Guinet C, Boutigny A-L, Fabre B, Frey P and Halkett F (2018) Defense Compounds Rather Than Nutrient Availability Shape Aggressiveness Trait Variation Along a Leaf Maturity Gradient in a Biotrophic Plant Pathogen. *Front. Plant Sci.* 9:1396. doi: 10.3389/fpls.2018.01396

Foliar pathogens face heterogeneous environments depending on the maturity of leaves they interact with. In particular, nutrient availability as well as defense levels may vary significantly, with opposing effects on the success of infection. The present study tested which of these factors have a dominant effect on the pathogen's development. Poplar leaf disks of eight maturity levels were inoculated with the poplar rust fungus *Melampsora larici-populina* using an innovative single-spore inoculation procedure. A set of quantitative fungal traits (infection efficiency, latent period, uredinia size, mycelium quantity, sporulation rate, sporulation capacity, and spore volume) was measured on each infected leaf disk. Uninfected parts of the leaves were analyzed for their nutrient (sugars, total C and N) and defense compounds (phenolics) content. We found that *M. larici-populina* is more aggressive on more mature leaves as indicated by wider uredinia and a higher sporulation rate. Other traits varied independently from each other without a consistent pattern. None of the pathogen traits correlated with leaf sugar, total C, or total N content. In contrast, phenolic contents (flavonols, hydroxycinnamic acid esters, and salicinoids) were negatively correlated with uredinia size and sporulation rate. The pathogen's fitness appeared to be more constrained by the constitutive plant defense level than limited by nutrient availability, as evident in the decrease in sporulation.

**Keywords:** constitutive defense, aggressiveness traits, plasticity, nutrient availability, rust fungus

## INTRODUCTION

Plant diseases and pests are responsible for crop damage that may account for up to 40% of yield losses worldwide (Boonekamp, 2012). Pathogen aggressiveness, the quantitative component of pathogenicity (Pariaud et al., 2009a), has been widely studied in order to forecast the evolutionary potential of plant pathogens and to better design disease control strategies (Lannou, 2012). Aggressiveness can be evaluated at different scales: by measuring the epidemic rate at the field scale or through elementary phenotypic traits at the individual scale (Pariaud et al., 2009a). The environment plays a key role in developmental processes of the pathogen and influences the individual's phenotype and the scale of fitness (Scheiner, 1993; Jarosz and Davelos, 1995).
Biotrophic fungal pathogens, such as rust and powdery mildew fungi, spend most of their life cycle within living plant tissues, so their fitness is directly affected by leaf physiology (Agrios, 2005; Gamm et al., 2011; Duplessis et al., 2013). On the one hand the leaf tissue provides the pathogen with the necessary nutrient resources for its development. On the other hand, it impedes pathogen development through a battery of plant defense compounds (Ullah et al., 2017).

Studies on interactions between plants and their foliar pathogens have shown that plant age, and especially leaf age, leads to different responses of the pathogen. Coleman (1986) postulated from a series of experiments that leaves pass through a peak of susceptibility to disease and pest attack. He proposed that for many stresses, including infection by biotrophic fungi, this peak coincides with changes in leaf properties during the transition from metabolite sink to source. This stage corresponds to a maximal nutrient availability. After this sink-source transition, free sugar content and thus nutrient resources decrease with increasing leaf maturity (physiological stage at a relative leaf age). Many studies support this view, demonstrating a maximum disease level before leaves are fully expanded and a decrease of disease severity on older leaves, which is also referred to as ontogenetic resistance (Ficke et al., 2002; Xavier et al., 2015; Farber and Mundt, 2017). In agreement, Merry et al. (2013) demonstrated that the peak in the severity of powdery mildew along a primary shoot of *Vitis vinifera* corresponds to the position of the first leaf that does not import carbohydrates. However, some studies have reported a contrasting pattern, e.g., Knott and Mundt (1991) found that infection efficiency of the wheat leaf rust fungus *Puccinia recondita* f. sp. *tritici* is 25% higher on older than on younger leaves in five wheat cultivars. The opposite trend may involve the dark side of plant-pathogen interaction: defense compounds. Their actions on a pathogen's development have been studied in several pathosystems, including rust fungi infecting poplar (Johnson and Kim, 2005), black currant (Vagiri et al., 2017), and willow (Hakulinen et al., 1999). Phenolic compounds present in leaves appear to highly influence the general plant resistance to fungi (Adandonon et al., 2017; Ullah et al., 2017). Most studies have focused on plant susceptibility to the disease (i.e., infection efficiency of the pathogen) or focused on overall disease severity. To dissect which fungal function is most influenced by leaf maturity and how, we propose to measure a set of elementary traits underlying components of aggressiveness.

Foliar pathogens of woody plants are good biological models to study the influence of leaf age on disease development, because they can infect leaves at various levels of maturity on the same host plant (Cellerino et al., 1978; Sharma et al., 1980; Turechek and Stevenson, 1998). Poplar is particularly well suited for studying the effect of leaf age on pathogens, since poplar plants grown from cuttings in a glasshouse produce new leaves at a constant rate, giving access to a wide range of leaves of increasing maturity on the same plant (Larson and Isebrands, 1971). These authors introduced the concept of leaf plastochron index (LPI) to refer to a position of a leaf along a poplar shoot. The youngest fully expanded leaf corresponds to LPI 1. Older leaves are numbered with increasing LPIs, according to their position along

the shoot. In most experimental conditions the sink to source transition in poplar is established around LPI 6 (Coleman, 1986), which thus defines the first mature leaf. The effect of leaf maturity on the development of poplar pathogens (*Marssonina brunnea*, *Melampsora medusae*, and *M. larici-populina*) has been widely studied (Cellerino et al., 1978; Sharma et al., 1980; Coleman, 1986; Giorcelli et al., 1996). In most of these experiments, leaves were spray-inoculated with spore suspensions. In these conditions, hundreds of spores compete for infection, and as a result, inoculum density has been shown to influence aggressiveness traits (Giorcelli et al., 1996; Pei et al., 2003; Pariaud et al., 2009a). A way to eliminate any possible interaction among spores and to remove the effect of lesion density on trait variation is to use a single-spore inoculation protocol. Moreover, the single-spore inoculation enables the examination of phenotypic variation at the individual lesion level and not by means of average trait values measured on a cohort of lesions.

The aim of this study was to assess the effect of leaf age, and the associated content in nutrients and defense compounds, on the variation in aggressiveness traits in the poplar rust fungus *M. larici-populina*. We quantified trait variation along the leaf maturity gradient using the same fungal genotype to establish a so-called reaction norm (Gavrilets and Scheiner, 1993). We measured the effect of LPI using a single-spore inoculation method for a suite of aggressiveness and morphological traits: infection efficiency, latent period, uredinia size, sporulation rate, sporulation capacity, in planta mycelium quantity, and urediniospore volume. In a first step, we examined the differences in pathogen development along a leaf maturity gradient in order to study the variation in fungal traits and to assess the correlations between them. In a second step, we estimated the influence of the most important leaf compounds (sugars, total carbon and nitrogen, and phenolics) on the variation in aggressiveness traits.

# MATERIALS AND METHODS

## Plant and Fungal Material, and Experimental Design

Characterization of quantitative traits was performed on excised leaf disks of *Populus deltoides* × *P. nigra* 'Robusta'. Poplar plants were grown from dormant cuttings in 10-L pots containing a sand-peat (1:1, v/v) mixture, with an initial fertilization of 3.5 g L⁻¹ CaCO₃ and 6 g L⁻¹ of slow release 13-13-13 N-P-K fertilizer (Nutricote® T100, Fertil). The plants were grown in a growth chamber regulated at 20/22°C (night/day temperatures) and 16 h photoperiod in order to ensure optimal growth, and were watered daily with deionized water. After four months, plants were about 1.2 m high and exhibited 25 to 30 fully expanded leaves.

On the same day, we harvested eight leaves of different maturity (LPI 6, 8, 10, 12, 14, 16, 18, and 20) from each of three plants (total of 24 leaves). From each poplar leaf, 24 disks (12-mm-diameter) were excised, and placed in flotation, each in a different cell culture plate. Each plate consists in 24 wells and thus contained one sample (leaf disk) of all leaves.
A plate thus corresponded to a complete statistical block: Twenty-four replicate plates were prepared, each containing the 24 combinations of LPI and plants. Within plates, the positions of the samples were randomized (Supplementary Figure S1). After disk excision, the remaining part of each leaf (about half a leaf) was stored individually in an aluminum wrapper at -80°C for sugar, nitrogen, carbon, and phenolic compound quantification. We hereafter refer to this plant material as "uninfected leaves."

Because the single-spore inoculation protocol is highly time-consuming, we were not able to inoculate more than six plates per day. Inoculations of the 24 plates were thus distributed over four consecutive days (each inoculation day formed a series of six plates). In other words, all leaf disks were excised the same day (day 1), but inoculated after 0 to 3 days depending on the series they belonged to.

Our experimental design (randomized complete block design) thus consisted in three experimental factors: LPI, plant and day of inoculation. Given an infection efficiency of approximately 25%, we expected six replicates per combination.

## Single-Spore Inoculation Protocol

Single-spore inoculation was performed with pure *M. larici-populina* isolate 98AG31. This isolate was collected in 1998 in Moÿ-de-l'Aisne (France) on *Populus trichocarpa* × *Populus deltoides* 'Beaupré' leaves and serves since as a reference isolate (Duplessis et al., 2011). Urediniospores were cryopreserved at -80°C at Institut National de la Recherche Agronomique Nancy (France) and multiplied on fresh poplar leaf disks before inoculation to maximize their germination and infection rates (Pei et al., 2002). To this aim, 1 mg of urediniospores was dispersed in 200 µL of water-agar (0.1 g L⁻¹). The resulting urediniospore suspension was applied as 2 µL droplets on the abaxial surface of two 30-mm-diameter leaf disks of 'Robusta'. Leaf disks were incubated by floating on deionized water in six-well polystyrene cell culture plates, with the abaxial surface upside, at 19 ± 1°C, and under continuous fluorescent light (25 µmol m⁻² s⁻¹). After a 10 to 13 days incubation period, the sporulating leaf disks were gently tapped over a microscope slide to release and collect urediniospores.

A germination test was performed just before inoculation to ensure the quality of urediniospores used. To this aim, a few hundred urediniospores were dispersed on the surface of a Petri dish containing agar (20 g L⁻¹). After an overnight incubation at 19 ± 1°C, the ratio of germinated/total urediniospores was evaluated under a light microscope (100× magnification) to determine the germination rate. This test ensured that nearly all spores used for the experiment were physiologically able to germinate (mean germination rate of 93%).

For single-spore inoculation leaf disks were inoculated one by one. One single urediniospore was picked with a human eyelash under a stereomicroscope (63× magnification) and deposited into a 5 µL water-agar droplet placed at the center of the

poplar leaf disk (abaxial surface). After inoculation the plates were incubated for 13 days in the same conditions as described above.

## Quantitative Trait Measurements

Several aggressiveness traits classically studied in plant-pathogen interactions were measured: infection efficiency, latent period, uredinia size, sporulation rate, and sporulation capacity (Lannou, 2012). In addition, we measured *in planta* mycelium quantity, which is less commonly studied.

* *Infection efficiency* is commonly defined as the probability that a urediniospore deposited on a receptive host surface produces a lesion (Lannou, 2012).

* *Latent period* is defined as the time interval between infection and the onset of sporulation from that infection (Lannou, 2012).

* *Uredinia size* is defined as the surface area that produces the urediniospores (Kolmer and Leonard, 1986; Robert et al., 2004) and has been suggested to be an indication of the pathogen capacity for host tissue colonization (Pariaud et al., 2009b; Lannou, 2012). However, the uredinia size is not a proxy of tissue colonization by *M. larici-populina*, as no necrosis is observed. Unlike necrotrophic pathogens, the only visual symptom is the opening in the leaf epidermis, through which urediniospores are extruded.

* *Sporulation rate* is defined as the number of urediniospores produced by the uredinia per unit of time (Leonard, 1969; Clifford and Clothier, 1974; Kardin and Groth, 1989): $\text{sporulation rate} = \text{urediniospore number}/(13 - \text{latent period})$.

* *Sporulation capacity* is defined as the amount of spores produced per unit area of sporulating tissue per unit of time (Robert et al., 2002; Dowkiw et al., 2003). It characterizes the flux of spores extruded by the uredinia and is calculated as the ratio between sporulation rate and uredinia size.

* *Mycelium quantity* which gives a proxy of the fungal colonization *in planta* is assessed using qPCR analysis on leaf disks (see below).

Beyond these six quantitative aggressiveness traits, we also measured a morphological trait:

* *Volume of urediniospores* was computed from length and width of urediniospores obtained from image analysis (see below).

We monitored the infections for 13 days, allowing measurements of latent period and infection efficiency. Emergence of uredinia was scored twice a day, at 9 a.m. and at 4 p.m. Infection efficiency was assessed for each leaf disk and was recorded as a binary variable: 1 if the deposited urediniospore produced one uredinia, 0 if not.

At day 13 (end of the monitoring period), leaf disks were harvested. We harvested all infected leaf disks (total: 158 samples out of 576 leaf disks inoculated, mean infection efficiency of 27.4%), giving from 4 to 12 biological replicates per LPI and per plant (mean: 6.58 replicates). The first step of the harvest consisted in separating the urediniospores produced
from the leaf disk. Infected leaf disks were placed in 2 mL Eppendorf tubes with 2 mL of Isoton® II isotonic buffered diluent (Beckman Coulter). The tubes were vigorously shaken for 20 s at 4 m s⁻¹ on a MP FastPrep®-24 homogenizer, in order to release all urediniospores from the uredinia. The second step consisted in measuring the size of the uredinia (or lesion size): photographs of the spore-free uredinia were taken under a stereomicroscope (25× magnification) coupled to a digital camera. Pictures were analyzed using ImageJ version 1.5i with a dedicated script (Supplementary Note S1). To measure the number of yellow pixels corresponding to the sporulating area, the Hue Saturation and Brightness was used: pixels with a hue under 37° were counted. Last, infected leaf disks were stored at -80°C and lyophilized before DNA extraction.

Measurements of both sporulation rate and urediniospore morphological parameters were performed on the urediniospore suspensions using an Occhio® Flowcell FC200+ optical morphogranulometer, which allows simultaneous particle counting and image analysis of the particles. Before counting, urediniospore suspensions were vortexed to disperse clusters of urediniospores. For each sample, counting was performed using the following settings: 0.15 mL of priming, 0.85 mL of volume analysis, and 7% of volume sampling (i.e., 7% within 0.85 mL were really analyzed, evenly distributed along the analyzed volume). Two filters parameter values were applied:

1. For urediniospore counting: bluntness between 0.6 and 1; length of ellipse minor axis between 20 and 60 µm; ellipse elongation between 0.18 and 0.79; and ISO eccentricity higher than 0.15,

2. For urediniospore size measurement: bluntness between 0.6 and 1; inertia Feret elongation between 0.25 and 1; ISO eccentricity higher than 0.1; minimum Feret diameter between 13.5 and 70 µm; and Feret diameter length between 23 and 200 µm.

These filters were applied on the raw data to compute sporulation rate expressed in spores per day, length, and width of urediniospores produced by each uredinia. The volume of urediniospores was computed following the formula of the volume of an ellipsoid (Philibert et al., 2011):

$$ \text{Spore volume} = \frac{4}{3}\pi \times \frac{\text{length}}{2} \times \left(\frac{\text{width}}{2}\right)^2 $$

## DNA Extraction and qPCR Quantitation of in planta Mycelium

DNA was extracted from infected poplar leaf disks using the DNeasy96® DNA plant kit (Qiagen). The Frozen Plant Tissue protocol (DNeasy 96 Plant Handbook, June 2015) was followed except that 20 1-mm-diameter glass beads were added to one 3-mm-diameter tungsten carbide bead in the 2 mL Eppendorf tubes. The 158 samples were randomly processed in two DNA extraction plates, which were stored at 4°C until qPCR analysis.

In order to correct for amplification differences among qPCR, three technical replicates were performed.

In order to quantify *M. larici-populina* mycelium in planta, quantitative Taqman® PCR was performed using new primers and probe designed for the rDNA ITS (Internal Transcribed Spacer). These primers and probe allowed a significant gain in sensitivity compared to a previously used primer pair, which was designed for a single-copy gene (Guinet et al., 2016). Reagent mix contained 3.1 µL of pure water, 6.5 µL of Mastermix®, 0.13 µL of primer ITS-Mlp-F (Primer sequence 5'-3': TGACTCTTTGTATAAACCATTACCC), 0.13 µL of primer ITS-Mlp-R (TCAAAGTTGCCTTTGAGATACG), and 0.13 µL of probe ITS-Mlp-P (6-FAM-TGCATTGTGGCCCGTCAAAA-BHQ1) per sample. A volume of 3 µL of extracted DNA was added to 10 µL of reagent mix. Standard samples consisted of calibrated quantities of plasmid DNA containing the cloned ITS amplicon. Two series of eight 10-fold dilutions of plasmid DNA (from 1.47 × 10⁷ to 1.47 fg µL⁻¹) were included in each qPCR run to set a standard curve. Quantitative PCR runs were performed on a QuantStudio 6 Flex (Life Technologies) machine and consisted in 10 min at 95°C followed by 45 cycles of 15 s at 95°C and 45 s at 62°C. DNA quantity was analyzed using the QuantStudio™ Real-time PCR software (version 1.1).

## Leaf Content Measurements: Total Nitrogen and Total Carbon Contents, Sugars, Phenolic Compounds

All leaf content measurements were conducted on the uninfected leaves, stored at -80°C just after disk excision.

Total nitrogen and total carbon concentrations in uninfected leaves were analyzed using an elemental analyser (Thermo Quest, Type NCS 2500). Measurements, data analysis, data record, and weighing record (Mettler Tolédo MT5 microbalance) were performed using the EAGER 200 software.

Soluble sugars were extracted from 20 mg of dry leaf powder with 1 mL of 70% (v/v) methanol/water for 10 min and then centrifuged at 17,000 g for 5 min at 4°C. This step was repeated twice and the resulting supernatants were pooled together. The supernatants, containing soluble sugars, were dried overnight with a vacuum evaporator (Maxi-Dry Plus; Hetomodel DW1, 0-110, Heto-HOLTEN A/S) to eliminate methanol. Dried extracts were rehydrated in 1.5 mL of distilled water, dissolved by sonication and then filtered at 0.2 µm (Acrodisc® Supor® filter; 0.2 µm, Van Waters Rogers). Undiluted aliquots of 20 µL were injected into a high-pressure liquid chromatography system (Dionex™ ICS-5000+ HPLC™) to determine soluble sugar composition. Soluble sugar separation was achieved on a Dionex™ CarboPac™ SA10 IC separation column (6 µm 250 × 4 mm; Thermo Scientific™). The flow rate was maintained at 1 mL min⁻¹, and the column temperature at 40°C. Individual carbohydrates were eluted from 0 to 30 min after injection. The sugar peaks were detected by light scattering (Sedex 45 ELSD system, Seder) then identified according to retention time and commercial standards. Sucrose, glucose, and fructose concentrations (% of dry weight) were determined by peak height and standard calibrations.

Phenolics were extracted from 30 mg dry leaf powder using the method described by Royer et al. (2013). Total phenolic content
was determined according to the Folin–Ciocalteu method. Briefly, each extract was diluted 1:16 in distilled water, and 20 µL of each diluted extract was distributed in a 96-well plate. Then 100 µL of Folin–Ciocalteu’s reagent (Sigma), diluted 1:10 in distilled water, and 80 µL of 7.5% sodium carbonate (w/v) were added. The microplate was incubated 30 min at 25°C and shaken 5 s every 10 min. Absorbance of the sample was measured with a microplate-reader (SynergyHT, BioTek) at 760 nm. Quantification of total phenolics was performed according to a calibration curve of gallic acid, and the results were expressed as milligram of gallic acid equivalent per gram of leaf dry weight.

In addition to total phenolic content, the concentration of major phenolics was analyzed on ultra-high performance liquid chromatography (U-HPLC) system (Shimadzu) equipped with a photo diode array detector (DAD) and a mass spectrometer. Each sample (3 µl) was separated on a C18 kinetex (100 mm × 2.1 mm) column (Phenomenex). The mobile phase consisted in 0.1% formic acid in ultra-pure water (solvent A) and 0.1% formic acid in methanol (solvent B). The molecules were eluted with a flow rate of 300 µL min⁻<fnref n="1" />, through a gradient elution from 1 to 50% B for 10 min, then to 99% B in 3 min, which was maintained for three additional minutes. The column was then re-equilibrated to 1% B for 4 min prior to the next run. Mass spectrometry analysis was carried out in ESI positive and negative modes. Quantification was performed by measuring the area under each peak at 280, 320, or 350 nm, depending on the lambda max of each molecule. Due to co-eluting signal at 280 nm, catechin quantification was realized by measuring the area under peak at *m/z* 289, corresponding to the [M-H]⁻ of the molecule. Experimental exact masses and MS fragments were compared to metabolomics databases (Respect<fnref n="1" />; MassBank<fnref n="2" />, DNP<fnref n="3" />) and data available in the literature in order to identify the nature of the metabolites. Sixteen phenolic compounds were identified. According to their structure, these 16 metabolites could be grouped into four phenolic subclasses, which are Flavan-3-ols (catechin, proanthocyanidin dimer), hydroxycinnamic acid esters (chlorogenic acid, caffeoyl shikimate, caffeoyl glucose isomers (×2), coumaroyl glucose isomers (×2)), flavonols (rutin, quercitrin), and salicinoids (salicine, salicortin, acetylsalicortin, homaloside D, tremulacin, and a putative populoside) (Supplementary Figure S2).

## Data Analysis

Principal component analyses (PCA) and linear models were used to explore correlations between phenotypic traits. Each independent phenotypic trait was analyzed separately using a generalized linear model evaluating effects of LPI, plant, day of inoculation, and the interaction between LPI and plant. To normalize the error, infection efficiency was analyzed with a binomial error family; latent period and sporulation capacity with a Gamma error family; uredinia size, spore volume, sporulation rate, and mycelium quantity with a Gaussian error family; and uredinia size and sporulation rate were square root

transformed. An analysis of variance (ANOVA) with an alpha error <5% was performed on each model. Differences of means were checked with Tukey tests. PCA on leaf characteristics was performed to explore their variations. Spearman correlation coefficients between phenotypic traits and leaf content were computed on mean phenotypic traits per LPI and plant based on 4 to 12 replicates. Spearman correlation coefficients' associated *p*-values were adjusted with a False-Discovery-Rate correction for multiple tests. All data analyses were performed using R version 3.2.5 (see details in the RMARKDOWN document, Supplementary Note S2).

## RESULTS

### Variation in *M. larici-populina* Quantitative Traits

To describe the variation in quantitative traits of *M. larici-populina* in this experiment, a PCA was performed on individual trait measurements assessed for each infected leaf disk (latent period, uredinia size, sporulation rate, sporulation capacity, mycelium quantity, and mean volume of produced spores). Infection efficiency cannot be included in the PCA analysis because it is set to one for infected leaf disk. Experimental factors (LPI, plant harvested, and day of inoculation) were considered as supplementary variables. The first principal component (PC1) explained 31% of the total inertia and was related to uredinia size, sporulation rate and to a lesser extent to sporulation capacity (Figure 1). PC2 and PC3 accounted for almost 20% and 17% of the total variance, respectively. PC2 was related to latent period which varied independently from uredinia size and sporulation rate, and PC3 was related to spore volume. However, there was large dispersal of individual measures that remained unexplained by experimental factors (Supplementary Note S2).

In accordance with PCA analysis, LPI was the main structuring variable in GLM analyses and had a significant effect on uredinia size ($P_{adjusted} = 1.22 \times 10^{-21}$), sporulation rate ($P_{adjusted} = 7.67 \times 10^{-15}$), and sporulation capacity ($P_{adjusted} = 1.40 \times 10^{-05}$) (Table 1). Uredinia size and sporulation rate increased with the LPI, while sporulation capacity decreased (Figure 2). Uredinia size exhibited a threefold increase (from 0.272 to 0.737 mm<fnref n="2" />) and sporulation rate a twofold increase (from 333 to 611 urediniospores day⁻<fnref n="1" />) from LPI 6 to 20 (Figure 2). Sporulation capacity was reduced more than twofold from LPI 6 to 20 (from 12,420 to 5,800 urediniospores day⁻<fnref n="1" /> mm⁻<fnref n="2" />) (Figure 2). The sporulation rate was also significantly, although less strongly, explained by the plant used for the experiment ($P_{adjusted} = 4.63 \times 10^{-05}$). Sporulation rate was significantly higher on plant No. 3 (mean rate = 526 urediniospores day⁻<fnref n="1" />) and lower on plant No. 2 (mean rate = 396 urediniospores day⁻<fnref n="1" />) (Figure 2), but no interaction was found between plant and LPI (Table 1). Other fungal traits varied but not in accordance with LPI, plant and day (Supplementary Figure S3).

We found no significant effect of the day of inoculation on any aggressiveness traits. As there was no difference in fungal traits across days of inoculation, we assumed that the incubation

<footnote n="1">¹http://spectra.psc.riken.jp/</footnote>

<footnote n="2">²http://massbank.jp/</footnote>

<footnote n="3">³http://dnp.chemnetbase.com/</footnote>
<figure><img src="image_22.png" /><figcaption>**FIGURE 1** | Principal component analysis of fungal phenotypic traits (latent period, uredinia size, sporulation rate, sporulation capacity, spore volume, and mycelium quantity, measured on 4 to 12 replicates per LPI and per plant); LPI, plant and day are supplementary variables. Colors correspond to variables' contribution to the first and second axes.</figcaption></figure>

<table><caption>**TABLE 1** | Summary of the analysis of variance of each experimental variable in generalized linear models (Chi-square test probability adjusted with a False-Discovery-Rate correction for multiple tests).</caption><thead><tr><th>Traits</th><th>LPI</th><th>Plant</th><th>LPI:Plant</th><th>Day</th></tr></thead><tbody><tr><td>Infection efficiency</td><td>0.459</td><td>0.661</td><td>0.346</td><td>0.165</td></tr><tr><td>Latent period</td><td>0.720</td><td>0.165</td><td>0.720</td><td>0.661</td></tr><tr><td>Uredinia size</td><td><strong>1.22 x 10<sup>-21</sup></strong></td><td>0.092</td><td>0.458</td><td>0.346</td></tr><tr><td>Spore volume</td><td>0.980</td><td>0.362</td><td>0.720</td><td>0.720</td></tr><tr><td>Sporulation rate</td><td><strong>7.67 x 10<sup>-15</sup></strong></td><td><strong>4.63 x 10<sup>-05</sup></strong></td><td>0.497</td><td>0.090</td></tr><tr><td>Sporulation capacity</td><td><strong>1.40 x 10<sup>-05</sup></strong></td><td>0.379</td><td>0.720</td><td>0.092</td></tr><tr><td>Mycelium quantity</td><td>0.866</td><td>0.980</td><td>0.996</td><td>0.720</td></tr></tbody></table>

P-values < 0.05 are indicated in bold characters.

conditions during the pre-inoculation period described above kept the leaf disks in their initial physiological state.

## Variation in Leaf Content Characteristics

To assess leaf content variations among the range of leaf maturity, a PCA was performed on leaf characteristics (Figure 3). Leaf characteristics were measured on one uninfected half-leaf of each plant and each LPI, so only plant and LPI supplementary variables were taken into account. PC1 explained almost 53%

of the total inertia and was highly related to sugars (sucrose, fructose, glucose) and phenolic compounds (total phenolics, and all the phenolic subclasses) variations. Explaining 18% and 12% of the inertia, PC2 and PC3 were associated to sucrose and nitrogen, respectively. The total carbon content of leaves did not vary among the range of leaf maturity studied. Leaf maturity, characterized by the content in phenolic compounds, appeared to structure more leaf characteristics than plant did.

## Correlation of Fungal Phenotypic Trait Variation and Leaf Content Characteristics

To assess how fungal and leaf data were correlated, pairwise Spearman's correlation coefficients and their adjusted *p*-values were calculated.

Spearman's coefficient for uredinia size and both sporulation rate and sporulation capacity were highly significant (rhô = 0.87 and -0.73, P adjusted = 3.37 x 10⁻⁶ and 4.66 x 10⁻⁴, respectively) (Figure 4), which is consistent with PCA results on individual traits. No other Spearman's coefficient appeared significant between fungal traits.

Spearman coefficients showed a negative correlation between total phenolic content and uredinia size and sporulation rate
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

**FIGURE 2** | Barplots of means and standard error calculated on three plants and 15 to 24 replicates per LPI (italic numbers at the bottom of the bars) of (A) uredinia size, (B) sporulation capacity, and (C) sporulation rate for each LPI; and (D) sporulation rate for each plant. Letters correspond to Tukey test results.

(rhô = -0.77 and -0.61, P adjusted = 1.62 × 10⁻⁴ and 1.07 × 10⁻², respectively), as well as a negative correlation between hydroxycinnamic acid esters content and uredinia size and sporulation rate (rhô = -0.83 and -0.68, P adjusted = 2.25 × 10⁻⁵ and 1.85 × 10⁻³, respectively), and a positive correlation between hydroxycinnamic acid esters content and sporulation capacity (rhô = 0.68, P adjusted = 1.85 × 10⁻³) (Figures 4, 5). Three of the four phenolic subclasses (salicinoids, hydroxycinnamic acid esters, and flavonols) showed the same trends as total phenolic content, with the exception that the correlation between sporulation rate and sporulation capacity, on one side, and salicinoids and flavonols, on the other side, were not significant (P adjusted > 0.05). In accordance with the PCA on leaf characteristics, all sugar contents were correlated with each other.

However, no correlation was found between fungal phenotypic traits and sugars, carbon and nitrogen contents. Salicinoids, flavonols, and hydroxycinnamic acid esters were all positively correlated to each other. Flavan-3-ols were, however, negatively correlated to the three other phenolic subclasses, although this was only significant with flavonols.

## DISCUSSION

In this study, we analyzed the variation in seven aggressiveness and morphological traits of the poplar rust fungus when inoculated on poplar leaves of varying maturity levels. A single-spore inoculation technique was used, which enabled
<figure><img src="image_27.png" /><figcaption>**FIGURE 3** | Principal component analysis on leaf chemical content (carbon, nitrogen, glucose, fructose, sucrose, total phenolics, flavan-3-ols, hydroxycinnamic acid esters, flavonols, and salicinoids), LPI and plant are supplementary variables. Colors correspond to variables' contribution to the first and second axes.</figcaption></figure>

us to study trait variance and covariance at the individual lesion level. Applying this method, we aimed at assessing the effect of leaf maturity on aggressiveness traits of *M. larici-populina*. In a second step, we explored correlations between the variation in fungal phenotypic traits and physiological features of the leaves in order to determine whether variation in phenotypic traits was influenced by nutritive or defensive compounds of leaves.

## Effect of Leaf Maturity on Aggressiveness and Morphological Traits

Regarding the variation in aggressiveness traits, our first result was a positive correlation between uredinia size and sporulation rate, whereby both increased with leaf maturity. That is, on more mature leaves, uredinia were wider with a higher production of urediniospores. Some authors found that the breaking down of sporulation rate into uredinia size and sporulation capacity enabled to explain differences in aggressiveness levels: an increase in sporulation rate can either be due to wider uredinia or higher sporulation capacity (flux of urediniospores) (Pariaud et al., 2009b). In our study, we found that the increase in sporulation rate is linked to wider uredinia, but not to higher sporulation

capacity. The correlation between uredinia size and sporulation rate has already been identified by several authors, including Dowkiw et al. (2003) on poplar rust, Pei et al. (2002) on willow rust and Robert et al. (2002) on wheat leaf rust. Based on our data, we add to this that sporulation capacity is not necessarily even as assumed by Robert et al. (2002). Indeed, sporulation capacity significantly decreased with leaf maturity, in contrast to uredinia size and sporulation rate. Even if the flux of urediniospores was lower on mature leaves, the amount of urediniospores produced was higher, leading to an increased aggressiveness on mature leaves compared to younger ones.

In this study, we also quantified individual variation in mycelium quantity through qPCR analysis. This trait has received little attention as a component of aggressiveness, mostly for methodological reasons. However, it can serve to measure the relative investment in sporulation vs. mycelium growth, the two main functions to which host extracted resources are allocated (Gilchrist et al., 2006; Yegorov et al., 2017). Notably, the quantity of *in planta* mycelium did not vary along the leaf maturity gradient and was correlated to neither sporulation rate nor lesion size. Hence, lesion size does not equate to mycelium
<figure><img src="image_28.png" /><figcaption>**FIGURE 4** | Pairwise Spearman correlation coefficients between leaf chemical content and mean fungal phenotypic traits per LPI and plant based on 4 to 12 replicates. Colors correspond to the sign and dot size to the strength of the correlation. Stars indicate significant correlations according to Student t-test adjusted with a False-Discovery-Rate correction for multiple tests.</figcaption></figure>

quantity. This makes sense as for biotrophic fungi like rusts, lesion (uredinia) size is only the size of the aperture through which urediniospores are extruded, and may thus not be directly linked to the size of the area within a leaf that is colonized by the mycelium.

Neither leaf maturity, nor the two other experimental factors (plant and day of inoculation) had an effect on latent period, infection efficiency, mycelium quantity and spore volume, even though substantial variances were observed for those traits. The latent period and spore volume varied independently from each other and independently from the main structuring traits (uredinia size and sporulation rate). Variation in latent period and spore volume might have been affected by an unidentified environmental factor. Even though inoculated leaf disks were incubated under controlled conditions of temperature and light, micro-scale variation may have occurred for these crucial parameters influencing fungal development and infection dynamics (Cochrane, 1959; Rotem et al., 1978; Trapero-Casas and Kaiser, 1992; Turechek and Stevenson, 1998; Turechek et al., 2001; Pap et al., 2013). Moreover, the initial position of the disk in the leaf might have influenced infection outcome

(Katsuya and Green, 1967; Johnston et al., 2016). Single-spore inoculation is more prone to inter-replicate variance, in contrast to spray-inoculation which buffers micro-environmental variations. Apart from this drawback, this method allowed us to accurately disentangle the variation in each trait and to properly assess putative correlations. In our study, the dissection of the infection process into elementary traits enabled us to assess that the leaf maturity gradient only influences the sporulation function.

## Disentangling Nutritive Limitation vs. Defense Reaction Effects

Two hypotheses were considered to explain variation in uredinia size and sporulation rate. On the one hand, sporulation could be enhanced by a higher nutrient availability in older leaves. On the other hand, higher levels of defense compounds could in contrast prevent fungal development in younger leaves.

Leaf sugars have long been known to be essential to the development of plant pathogenic fungi (Bilgrami and Verma, 1978; Madan and Thind, 1998). Using *Populus deltoides* leaves
<figure><img src="image_29.png" /><figcaption>**FIGURE 5** | Bivariate plots (LPI and plant mean ± standard error calculated on three plants and 15 to 24 replicates per LPI) showing the variation in (A) uredinia size and (B) sporulation rate according to the phenolic subclass hydroxycinnamic acid esters content. Colors stand for LPI and R2 gives the coefficients of determination.</figcaption></figure>