Population BiologyFree Access icon

The Melampsora americana Population on Salix purpurea in the Great Lakes Region Is Highly Diverse with a Contributory Influence of Clonality

    Affiliations
    Authors and Affiliations
    • Chase R. Crowell1
    • Dustin G. Wilkerson2
    • Mariami Beckauri1
    • Ali R. Cala1
    • Patrick W. McMullen1
    • Stephen Mondo3
    • William Andreopoulos4
    • Anna Lipzen3
    • Kathleen Lail3
    • Mi Yan3
    • Vivian Ng3
    • Igor V. Grigoriev3 5
    • Lawrence B. Smart2
    • Christine D. Smart1
    1. 1Plant Pathology & Plant-Microbe Biology Section, School of Integrative Plant Science, Cornell AgriTech, Cornell University, Geneva, NY 14456
    2. 2Horticulture Section, School of Integrative Plant Science, Cornell AgriTech, Cornell University, Geneva, NY 14456
    3. 3U.S. Department of Energy Joint Genome Institute, Walnut Creek, CA 94598
    4. 4Department of Computer Science, San Jose State University, San Jose, CA 95192
    5. 5Department of Plant and Microbial Biology, University of California, Berkeley, CA 94720

    Published Online:https://doi.org/10.1094/PHYTO-05-21-0201-R

    Abstract

    Shrub willows (Salix spp.) are emerging as a viable lignocellulosic, second-generation bioenergy crop with many growth characteristics favorable for marginal lands in New York State and surrounding areas. Willow rust, caused by members of the genus Melampsora, is the most limiting disease of shrub willow in this region and remains extremely understudied. In this study, genetic diversity, genetic structure, and pathogen clonality were examined in Melampsora americana over two growing seasons via genotyping-by-sequencing to identify single-nucleotide polymorphism markers. In conjunction with this project, a reference genome of rust isolate R15-033-03 was generated to aid in variant discovery. Sampling between years allowed regional and site-specific investigation into population dynamics, in the context of both wild and cultivated hosts within high-density plantings. This work revealed that this pathogen is largely panmictic over the sampled areas, with few sites showing moderate genetic differentiation. These data support the hypothesis of sexual recombination between growing seasons because no genotype persisted across the two years of sampling. Additionally, clonality was determined as a driver of pathogen populations within cultivated fields and single shrubs; however, there is also evidence of high genetic diversity of rust isolates in all settings. This work provides a framework for M. americana population structure in the Great Lakes region, providing crucial information that can aid in future resistance breeding efforts.

    Willow rust caused by Melampsora spp. is a devastating disease of shrub willow (Salix spp.) bioenergy production in the northeast United States (Karp et al. 2011; Smart and Cameron 2008). Heavy infection of these highly polycyclic pathogens results in premature defoliation, which can lead to 50% yield loss in susceptible hosts (Serapiglia et al. 2013; Verwijst 1990). The use of fungicides has been effective for willow rust disease control; however, they are expensive, and their use is not sustainable (McCracken and Dawson 1992, 2003). For these reasons, deployment of resistant Salix cultivars is essential for economic and environmental sustainability of this long-lived, perennial crop. Scientific efforts are accumulating to achieve this goal, including the development of genetic resources for the willow host Salix purpurea (Zhou et al. 2020) (https://phytozome.jgi.doe.gov/), as well as the discovery of quantitative trait loci in Salix species for resistance to willow rust in New York, Pennsylvania, and West Virginia (Carlson and Gouker et al. 2019). However, little is known about the Melampsora spp. pathogen populations. Investigations into genetic diversity can provide information about the genetic variation present in the sampled populations, enable hypotheses involving interspecific competition, and describe community structure (Hughes et al. 2008). In the context of plant–pathogen interactions, estimating genetic diversity may provide valuable insight into the expected success of deployed plant host resistance. This information would be very valuable in the S. purpureaMelampsora spp. pathosystem, because host resistance is the primary mechanism of disease control.

    Previous work has identified the species Melampsora americana as the primary contributor to disease on S. purpurea in New York (Crowell et al. 2020; Kenaley et al. 2014). M. americana is a macrocyclic and heteroecious rust alternating on the balsam fir (Fig. 1) (Kenaley et al. 2014). Little is known about the diversity of this species; however, inferences can be made by examining literature from related rust species. Investigations into population structure have been performed on willow rust species in other locations, including M. larici-epitea in Sweden (Samils et al. 2001), members of the M. larici-epitea species complex in Germany (Bubner et al. 2014), M. epitea in England (Pei et al. 2000), and various Melampsora species in western North America (Bennett et al. 2011). These studies explored genotypic diversity by using amplified fragment length polymorphism (AFLP) and rDNA phylogeny, which have limited resolution, and focused primarily on Melampsora spp.–Salix spp. host pathogenicity barriers. However, Samils et al. (2001) investigated the relative role of sexual and asexual reproduction of M. larici-epitea at three sites in Europe by using AFLP markers. In two sites in Sweden, they found high diversity of AFLP polymorphisms with no evidence for population structure. However, the third site in Ireland showed less diversity of polymorphism, indicating a moderate level of population differentiation. Although it is presumed that the sexual stage of M. americana is essential in the Great Lakes region of North America for completion of its life cycle, it has been noted that stem-infecting forms exist in Melampsora spp. rust infecting S. viminalis in Europe that forgo sexual reproduction on their alternate host (Pei et al. 1995, 2000). Determining the presence of sexual reproduction in comparison with clonal spread in a field or on a single shrub will further our understanding of the epidemiology of this pathogen.

    Fig. 1.

    Fig. 1. Life cycle diagram of Melampsora americana. Alternation of host is denoted by dotted lines, and approximate time of year is indicated in boxes.

    Download as PowerPoint

    Although little is known about population biology of M. americana, greater progress has been made in understanding M. larici-populina in Europe and Asia (Barrès et al. 2008; Gérard et al. 2006; Pei et al. 2007). These studies used random amplified polymorphic DNA markers, microsatellites, and AFLP markers to investigate the presence of genetic structure by location of collection. In all cases, frequent gene flow was observed, and there was limited evidence of geographic separation. However, Pei et al. (2007) found one site in northeast England that was genetically differentiated from the other populations. If the M. americana population shows similar patterns of genetic diversity to M. larici-populina in Europe, it is likely that the M. americana population is panmictic with frequent gene flow events; however, there may also exist locations with isolated populations.

    The goal of this project was to determine the genotypic diversity of the M. americana pathogen at local site and regional scales in the growing seasons of 2015 and 2017 collected from wild and cultivated S. purpurea (with some unknown host species). Although previous studies in willow and poplar rust have relied on AFLP, random amplified polymorphic DNA, and rDNA marker sets, reduced representation genome sequencing methods such as genotyping-by-sequencing (GBS) have since been developed and provide higher marker density and greater genetic resolution of individuals within a population (Aoun et al. 2020; Elshire et al. 2011). We hypothesize that the M. americana population in our sampled regions will show a high level of genetic diversity indicative of a frequently sexually reproducing population. Additionally, we hypothesize that there is some level of genetic structure by isolate collection site. This information will provide an understanding of pathogen genetic diversity at a regional and field level and will enable breeding programs to test for resistance against rust isolates representative of the breadth of pathogen diversity in this region.

    MATERIALS AND METHODS

    M. americana isolate collection and DNA extraction.

    A total of 182 M. americana isolates were collected from wild and cultivated S. purpurea in the summer and fall of 2015 and 2017. Ten symptomatic rust-infected leaves per shrub were collected in the field, stored in plastic bags in a cooler, and returned to the lab for single pustule isolations. Three uredospore pustules from each set of 10 leaves were used to generate three separate single-genotypic isolates by performing three successive rounds of single-pustule inoculations according to methods described by Crowell et al. (2020). Isolates (Supplementary Table S1) were maintained on S. purpurea ‘Fish Creek’ in water agar Petri dishes according to methods described by Crowell et al. (2020). Briefly, in 2015 isolates (n = 101) were collected from six states with a combination of wild and cultivated S. purpurea hosts (instances where host species were not certain are indicated in Supplementary Table S2). Rust isolates in 2017 (n = 81) were collected from cultivated S. purpurea at four fields (Supplementary Table S3). One field was a S. purpurea F2 mapping population in Geneva, New York (named Geneva NY 1) (Carlson and Gouker et al. 2019), and the remaining three were replicated plantings of a S. purpurea association mapping population in Geneva, New York (named Geneva NY 2), Portland, New York, and Morgantown, West Virginia (Carlson and Gouker et al. 2019). Sampling was opportunistic from wild Salix hosts and followed a W-sampling scheme within cultivated fields. All instruments used in isolation and collection of rust uredospores were sterilized between isolates with 70% ethanol. After single-pustule isolations, each individual isolate was grown on approximately 10 detached S. purpurea leaves, uredospores were collected with a cyclone spore collector, and DNA was extracted according to the protocol described by Crowell et al. (2020).

    Rust genome assembly.

    A M. americana genome assembly was generated to serve as a reference for single-nucleotide polymorphism (SNP) discovery. To obtain high-quality DNA for genome sequencing, M. americana isolate R15-033-03 was paintbrush inoculated onto detached leaves of S. purpurea ‘Fish Creek’ via methods described by Crowell et al. (2020). A total of 40 μg of high-molecular-weight genomic DNA (gDNA) was isolated via a protocol described by Duplessis et al. (2011) and Crowell et al. (2020) and submitted to the Joint Genome Institute (JGI) for whole genome sequencing. Ten micrograms of genomic DNA was sheared to fragments of approximately 30 kb with Megaruptor3 (Diagenode, Denville, NJ). The sheared DNA was treated with exonuclease to remove single-stranded ends and DNA damage repair mix, followed by end repair and ligation of blunt adapters with a SMRTbell Template Prep Kit 1.0 (Pacific Biosciences, Menlo Park, CA). The final library was size selected with BluePippin (Sage Science, Beverly, MA) at 15 kb cutoff size and purified with AMPure PB beads (Pacific Biosciences, Menlo Park, CA). PacBio Sequencing primers were then annealed to the SMRTbell template library, and sequencing polymerase was bound to them with a Sequel Binding kit 2.0 (Pacific Biosciences, Menlo Park, CA). The prepared SMRTbell template libraries were then sequenced on a Pacific Biosystems Sequel sequencer with v3 sequencing primer, 1 M v2 SMRT cells, and version 2.0 sequencing chemistry with 6- and 10-h movie run times. Genome assembly was generated in Falcon from the mitochondria-filtered PacBio reads and polished with Arrow SMRTLink version 5.0.1.9578 (https://github.com/PacificBiosciences/GenomicConsensus).

    Additionally, 10 μg of total RNA was extracted with a Spectrum Plant Total RNA Kit (MilliporeSigma, Burlington, MA) from diseased willow leaf tissue pooled from six growth stages (4, 5, 6, 7, 8, 9, and 10 days postinoculation). This was done to aid in annotation of the draft genome. Plate-based RNA sample prep was performed on a PerkinElmer Sciclone NGS robotic liquid handling system (PerkinElmer, Waltham, MA) with Illumina’s TruSeq Stranded messenger RNA HT sample prep kit via poly-A selection of messenger RNA according to the protocol outlined by Illumina (https://support.illumina.com/). Total RNA starting material was 1 μg per sample, and eight PCR cycles were used for library amplification. The prepared libraries were quantified with KAPA Biosystems’ next-generation sequencing library quantitative PCR kit and run on a Roche LightCycler 480 real-time PCR instrument (Roche, Basel, Switzerland). The quantified libraries were then prepared for sequencing on the Illumina HiSeq sequencing platform with a TruSeq paired-end cluster kit version 4. Sequencing of the flow cell was performed on the Illumina HiSeq 2500 sequencer with HiSeq TruSeq SBS sequencing kits version 4, according to a 2 × 150 indexed run recipe. Contigs <1,000 bp were excluded and RNA sequence mapping required removing S. purpurea plant from the RNA reads. Quality and artifact filtered reads were assembled into consensus sequences in Trinity version 2.3.2 (Grabherr et al. 2011). The genome was annotated on the JGI Annotation pipeline (Grigoriev et al. 2014).

    GBS analysis.

    Restriction enzyme optimization of GBS was required for M. americana before analysis of DNA from all isolates (Elshire et al. 2011). Rust isolate R15-033-03 (Crowell et al. 2020) was increased on ‘Fish Creek’ detached leaves, and DNA was extracted according to the protocol described by Crowell et al. (2020). The resulting 10 μg of pooled high-quality gDNA of this isolate was used for GBS restriction enzyme optimization by the University of Wisconsin-Madison Biotechnology Center with ApeKI.

    The M. americana rust isolates from 2015 and 2017 were propagated, and DNA was isolated according to the method of Crowell et al. (2020). Three plates containing 25 µl of ≥20 ng/µl gDNA per sample were submitted to University of Wisconsin-Madison Biotechnology Center for 96-plex GBS library preparation and next-generation sequencing. Two of these were sequenced via Illumina HiSeq2500 yielding single-end 100-bp reads containing isolates from 2015 and 2016, and the third plate was sequenced via NovaSeq6000 generating paired-end 150-bp reads containing isolates from 2017. This change in sequencing method reflected changes in the Biotechnology Center’s sequencing infrastructure that occurred between plate submissions.

    Barcode sequences were trimmed, and genotype calling was performed via the Tassel GBSv2 pipeline with default parameters (Bradbury et al. 2007; Glaubitz et al. 2014). Because sequencing products differed between the HiSeq2500 and NovaSeq6000, only the forward reads of the latter method were used for consistency. GBS reads were aligned to the M. americana R15-033-03 v1.0 reference genome (generated above) for variant discovery via the Burrows-Wheeler algorithm bwa-aln version 0.7.17 with default parameters (Li and Durbin 2010). Tassel 5 was used to remove individuals with >40% missing data. The following filtering steps on genotype sites were performed in Vcftools version 0.1.17 (Danecek et al. 2011): Genotypes with <5× coverage were set to missing, minor allele frequency cutoff was set to ≥0.05, mean site depth was between >10× and <50×, missing data cutoff was <20%, and indels were removed. This resulting vcf file was used to generate three separate datasets: isolates collected from 2015, isolates collected from 2017, and a pooled set of isolates from 2015 and 2017. These datasets were filtered again to remove sites with >40% missing data and minor allele frequency > 0.05 with Tassel 5 and VCFtools. All scripts and variant resources are available in the following GitHub repository (https://github.com/crc255/M_americana_population_biology_data).

    Population diversity analysis.

    Genotypic diversity analysis was conducted in R version 4.0.3 (R Core Team 2013) in RStudio version 1.4.1103 (RStudio Team 2020). Isolate a priori population designation for rust isolates collected in 2015 was determined with an 80-km geographic distance cutoff. The command ‘distm’ in the package ‘geosphere’ version 1.5-10 (https://CRAN.R-project.org/package=geosphere) was used to generate a geographic distance matrix between each isolate with latitude and longitude coordinates of collection sites. The command ‘hclust’ and ‘cutree’ in the package ‘stats’ version 2.15.3 (R Core Team 2013) was used for hierarchical cluster analysis based on the ‘distm’ geographic distance matrix with cutoff distance of 80 km for any two isolates to be considered in the same population. The 80-km cutoff was determined by visual best fit for the sampling sites assayed. A geographic map of isolates collected was generated via the ‘map’ function in ‘maps’ version 3.3.0 (https://CRAN.R-project.org/package=maps) with the geographic groups discovered based on the 80-km distance cutoff (Fig. 2).

    Fig. 2.

    Fig. 2. Population diversity of Melampsora americana isolates collected in 2015. A, Neighbor joining tree of 2015 isolates without clone correction. Each node represents a single isolate, and the axis indicates genetic distance between isolates, calculated via Provesti’s distance and 1,000 bootstrap replicates. B, Map of M. americana isolate collection sites. C, Plot of principal components 1 and 2 (PC1 and PC2) of clone corrected data. Percentage variation explained by each principal component is indicated on the axis. Size of point correlates with the number of isolates within a clonal lineage. Colors/shading across A, B, and C are consistent and indicate the geographic grouping from which the isolate was collected. The two largest clonal lineages are indicated by asterisks. Examples of clonal lineages with an isolate from a different geographic group from the others in that lineage are indicated by hashmarks.

    Download as PowerPoint

    Pairwise identity by state (IBS), defined as the proportion of alleles shared at nonmissing sites, was used with code provided by Gregory Vogel (Vogel et al. 2021) to approximate the sequencing error rate. A clonality cutoff was determined by comparing proportions of shared alleles between 11 technical replicates submitted for GBS. This approximation was used to establish a percentage identity cutoff for clonal lineages in the larger SNP dataset. Once this cutoff was determined, IBS was calculated across all individuals, and this cutoff was used to cluster isolates into clonal lineages. A clone corrected dataset was generated of one isolate per clonal lineage. We chose the isolate by computing the percentage missing data for each isolate and selecting the isolate with the least missing data per clonal lineage.

    Vcf files were read into R in the package ‘vcfR’ version 1.12.0 (Knaus and Grünwald 2017), and phylogenetic neighbor joining trees with Provesti’s distance with 1,000 bootstrap replicates were generated via the ‘aboot’ function in the package ‘poppr’ version 2.9.0 (Kamvar et al. 2014). Principal components of clone corrected datasets were calculated via the ‘glPca’ function in the package ‘adegenet’ version 2.1.3 (Jombart and Ahmed 2011). Determination of population separations was determined via Fst for isolates collected 2017 with a priori population designation (Geneva NY 1, Geneva NY 2, Portland NY, and Morgantown WV). Pairwise Fst values were calculated via the ‘stampfst’ function in the package ‘StAMPP’ version 1.6.1 (Pembleton et al. 2013). Pairwise Jost D was calculated for isolates collected in 2017 via the ‘pairwise_d’ function in the R package ‘mmod’ version 2.9.0 (Winter 2012). Sequencing data from all isolates collected in 2015 were subsampled to 39 samples from geographic groups containing isolates collected from both wild and cultivated hosts. Pairwise Fst between isolates collected from wild or cultivated hosts was calculated as described previously. The ‘poppr’ function in the R package ‘poppr’ version 2.9.0 (Kamvar et al. 2014) was used to generate summary genetic diversity statistics for rust isolates in collected in 2017. Discriminant analysis of principal components was performed via the ‘dapc’ function, and stacked compoplots were made via the function ‘compoplot’ in the package ‘adegenet’ version 2.1.3 (Jombart and Ahmed 2011).

    RESULTS

    Rust genome assembly and annotation.

    High-quality gDNA of M. americana isolate R15-033-03 was used for PacBio sequencing, genome assembly, and annotation at JGI. The total assembly size of the M. americana genome was 112.35 Mb over 389 scaffolds (Table 1). The average coverage depth was 128.49× with a maximum scaffold size of 3.25 Mbp. The N50 value was 34, and L50 value was 1.07 Mbp. The genome size is similar to the closely related poplar rust M. larici-populina genome (Duplessis et al. 2011); however, M. americana has significantly fewer repeats, with 27% of the genome compared with 44% in M. larici-populina (Table 1). The Core Eukaryotic Genes Mapping Approach (CEGMA) was used to assess completeness of the genome sequence. CEGMA capture was 97.16%. The genome was resequenced to improve capture; however, CEGMA values were unchanged. RNA sequencing of infected leaves was used to improve gene annotation, resulting in a total of 15,984 gene models. After contaminating host S. purpurea reads were removed, 378,939,217 expressed sequence tags were identified in 60,136 clusters after assembly, with 93% of reads aligned to the generated rust fungal genome. Genome annotation was completed, and the following summary information is available through the JGI MycoCosm platform: total carbohydrate-active enzymes, 408; total peptidases, 208; total secondary metabolism clusters, five (four nonribosomal peptide synthase or synthase-like, one polyketide synthase); total transcription factors, 456; and total transporter proteins, 810. This genome was made available through the JGI fungal genome portal MycoCosm as M. americana R15-033-03 v1.0 (https://mycocosm.jgi.doe.gov/Melame1/Melame1.home.html) and deposited at DDBJ/ENA/GenBank under accession number JAIFHY000000000.

    TABLE 1. Summary statistics of Melampsora americana reference genome R15-033-03 sequencing and RNA read mapping to genomea

    M. americana isolate GBS analysis.

    Rust uredospores of 182 M. americana isolates were successfully collected, maintained on ‘Fish Creek’, and used for DNA extraction. GBS data from all M. americana isolates processed through the Tassel 5 pipeline resulted in 177,177 variants (Bradbury et al. 2007; Glaubitz et al. 2014). After filtering, 5,865 variants were retained, resulting in a dataset containing high-quality SNPs shared between isolates collected from 2015 and 2017. This dataset was divided into two subdatasets with only 2015 or 2017 isolates, which were again filtered, resulting in 4,717 and 5,055 variants, respectively.

    Latitude and longitude data were used to determine geographic groupings based on location of collection of each isolate from 2015. The 10 geographic groups are summarized as follows (Fig. 2B): isolates collected in groups 1, 3, 4, 5, 7, and 8 were located in New York (n = 73); group 2 included all isolates from Pennsylvania (n = 6), group 6 included all isolates from Vermont (n = 15), group 9 isolates were from West Virginia (n = 5), and group 10 included all isolates from in Michigan (n = 2). Isolates collected in 2017 were grouped by the field of collection: Geneva, New York 1 (n = 15), Geneva, New York 2 (n = 21), Portland, New York (n = 20), and Morgantown, West Virginia (n = 25).

    A total of 11 technical replicates were included on the plates submitted for GBS to estimate error rate. Pairwise IBS was calculated for each technical replicate, resulting in an approximate overall error rate of 1.8%. The error rate for technical replicates of rust isolate R15-033-03 included on each of the three plates was 2.7%, with the greatest error rate between any two technical replicates of 3.4%. Therefore, a relaxed error rate cutoff of 4% was used to determine clonality across the full SNP dataset. One technical replicate was retained per genotype for the following analyses.

    Genetic diversity of isolates collected in 2015.

    Isolates collected in 2015 spanned the widest geographic range, including all 10 geographic groups described previously (Fig. 2B). A neighbor joining tree using Provesti’s genetic distance and 1,000 bootstrap replicates showed some clustering of isolates based on location of collection (Fig. 2A). However, no monophyletic clades contained all isolates from a single geographic region or group of locations. Isolates from 2015 also included M. americana collected from both wild (n = 23) and cultivated (n = 78) S. purpurea (Supplementary Table S2). These isolates were subsampled to geographic groups 3 and 4, where rust was collected from both wild (n = 19) and cultivated (n = 20) hosts. Pairwise Fst calculated on these subsampled data resulted in no population separation (Fst = 0.017) between wild and cultivated willow from these locations (Supplementary Table S4). Additionally, a neighbor joining tree using Provesti’s genetic distance showed no clustering of M. americana isolated from wild S. purpurea when compared with cultivated S. purpurea (Supplementary Fig. S1).

    Of the 101 total isolates from 2015, clone correction identified 74 distinct genotypes, of which 79% (59/74) were represented by a single isolate (Table 2). The largest clonal lineage contained nine isolates from five different S. purpurea shrubs from geographic group 1 and the second largest clonal group contained five isolates from three different S. purpurea shrubs from geographic group 9 (Fig. 2A). All other clonal lineages contained three or fewer isolates (Supplementary Table S2). We attempted three isolations of M. americana from each of 60 different shrubs, with a single successful isolation from 26 shrubs, two isolates from 27, and success of all three isolations from seven shrubs. Of the 27 shrubs yielding two isolates, the pairs of isolates from only six of those shrubs were clonal. The isolates from three of the seven shrubs yielding three isolates were clonal.

    TABLE 2. Summary data for clonality of Melampsora americana isolates within a single shrub for years 2015, 2017, and 2015 and 2017 together

    Principal component analysis of the clone corrected dataset separated three isolate genotypes from geographic group 8 from all other isolate genotypes, with principal component 1 explaining 3.18% of the variation (Fig. 2C). These three isolate genotypes (R15-072-01, R15-073-02, and R15-075-02) were collected in a cultivated shrub willow nursery on hosts 05-OSU-06, 04-BN-046, and 04-BN-04, respectively. Principal component 2 explained 3.11% of the variation and separated three isolate genotypes (R15-023-01, R15-063-03, and R15-020-02) from groups 3, 6, and 10, respectively, from the others. These isolates were collected from very distant locations on hosts 94001, 94006, and ‘Fish Creek’, respectively. Genotypes representing the two largest clonal isolate lineages do not separate from the other genotypes on PC1 or PC2. Principal component 3 explained 2.6% of the variation and separated the isolate genotype representing the second largest clonal lineage, collected in geographic group 9, from the remaining genotypes (Supplementary Fig. S2). Principal component 4 explained 2.4% of the variation and separated two isolate genotypes (R15-043-01 and R15-031-02) from groups 4 (on ‘Fish Creek’) and 3 (on an unknown wild Salix genotype), respectively.

    Genetic diversity of isolates from 2017.

    Isolates were collected from four S. purpurea breeding populations in 2017 (Fig. 3B). A neighbor joining tree identified isolate genotype clusters based on location of collection; however, many of these clades contains one or more isolates from different locations (Fig. 3A). This finding is corroborated by Fst calculations in which fields Geneva, New York 1, Geneva, New York 2, and Portland, New York showed low pairwise Fst values (0.042 to 0.059) and low Jost D values (0.020 to 0.028), suggesting no population separation by site (Table 3). However, field Morgantown, West Virginia shows some evidence of moderate population separation, with Fst values from 0.256 to 0.303 and Jost D values ranging from 0.112 to 0.137 for all pairwise calculations. Discriminant analysis of principal components (DAPC) was performed on these isolates to assess population structure without a priori designation, identifying six clusters explained by the data (Supplementary Fig. S4). The DAPC was used to explore genetic composition explained by these six clusters and was compared with genetic composition described with the a priori population designation based on field of collection (Supplementary Fig. S5).

    Fig. 3.

    Fig. 3. Population diversity of Melampsora americana isolates collected in 2017. A, Neighbor joining tree of 2017 isolates without clone correction. Each node represents a single isolate, and the axis indicates genetic distance between isolates, calculated via Provesti’s distance and 1,000 bootstrap replicates. B, Map of M. americana isolate collection sites. C, Plot of principal components 1 and 2 (PC1 and PC2) of clone corrected data. Percentage variation explained by each principal component is indicated on the axis. Size of points correlates with the number of isolates within a clonal lineage. Colors/shading across A, B, and C are consistent and indicate the geographic grouping from which the isolate was collected. The two largest clonal lineages are indicated by asterisks. Examples of clonal lineages with at least one isolate from a different geographic group from the others in that lineage are indicated by hashmarks.

    Download as PowerPoint

    TABLE 3. Pairwise Fst and pairwise Jost D calculations for Melampsora americana isolates collected in 2017 by field of collection

    Genetic diversity indices were calculated for all isolates collected in 2017 (Table 4). The isolate population collected from Morgantown, West Virginia showed significantly lower genotypic diversity than the other three populations according to the Shannon-Weiner diversity index (Shannon 2001), Stoddard and Taylor’s index (Stoddart and Taylor 1988), the corrected Simpson index (Simpson 1949), and Evenness (Grünwald et al. 2003; Ludwig and Reynolds 1988; Pielou 1975). Isolate populations from Geneva, New York 1, Geneva, New York 2, and Portland, New York showed similar levels of high genetic diversity across these statistics.

    TABLE 4. Summary statistics of genotypic diversity of Melampsora americana isolates collected in 2017 from select geographic groups

    Of the 81 isolates collected in 2017, there were only 36 unique genotypes, with 69% (25/36) of the genotypes represented by only a single isolate, in contrast to one clonal genotype represented by 23 of the isolates (Table 2). Among those 23 isolates, 22 were collected from 10 S. purpurea shrubs from Morgantown, West Virginia, and one isolate was collected from a S. purpurea shrub in Geneva, New York 2 (Fig. 3A). The second largest clonal lineage consisted of six isolates collected from four Portland, New York shrubs. Of the 16 clonal lineages containing more than one isolate, three included at least one isolate collected from a different field than the other isolates in that lineage (Supplementary Table S3). We attempted three isolations of M. americana from each of 39 different shrubs, with a single successful isolation from nine shrubs, two isolates from 16, and success of all three isolations from 13 shrubs. Of the 16 shrubs yielding two isolates, the pairs of isolates from only four of those shrubs were clonal. The isolates from 10 of the 13 shrubs yielding three isolates were clonal.

    Principal component analysis separated the genotype representing the largest clonal lineage, as well as two closely clustering genotypes, R15-025-02 from Geneva, New York 2 and R15-046-02 from Morgantown, West Virginia with principal component 1 explaining 6.9% of the variation (Fig. 3C). Principal component 2 explained 6.1% of the variation and separated those same three genotypes from most others. Principal component 3 explained 4.3% of the variation and showed a continuous spread of distribution (Supplementary Fig. S2). Principal component 4 explained 4.1% of the variation and separated two genotypes (R15-022-02 and R15-037-01) collected in Geneva, New York 2 and Portland, New York from the remaining genotypes.

    Combined genetic diversity analysis of isolates from 2015 and 2017.

    A neighbor joining tree was generated including all isolates from both 2015 and 2017 (Supplementary Fig. S3). Fields sampled in 2017 were grouped into the 10 locations described in the 2015 isolate data for comparison. Repeated isolation in successive years of collection was limited to groups 1, 8, and 9. Some clustering of isolates based on location of collection can be seen in this neighbor joining tree. Separately, a neighbor joining tree containing all isolates collected was generated instead, highlighting year of collection (Fig. 4). This tree did not show evidence of clustering of isolates based on year of collection.

    Fig. 4.

    Fig. 4. Neighbor joining tree based on Provesti’s distance and 1,000 bootstraps of Melampsora americana isolates collected in 2015 and 2017 without clone correction. Colors/shading indicate year of collection of isolate, and the two largest clonal lineages are indicated by asterisks.

    Download as PowerPoint

    Of the 182 total isolates collected in 2015 and 2017, there were 110 unique genotypes, with 76% (84/110) of these represented by a single isolate (Table 2). The largest clonal lineage from the combined dataset was the largest clonal lineage collected in 2017, containing 22 isolates from Morgantown, West Virginia and one from Geneva, New York 1 (Fig. 4). We attempted three isolations of M. americana from each of 99 different shrubs, with a single successful isolation from 35 shrubs, two isolates from 43, and success of all three isolations from 20 shrubs. Of the 43 shrubs yielding two isolates, the pairs of isolates from only 10 of those shrubs were clonal. The isolates from 13 of the 20 shrubs yielding three isolates were clonal.

    DISCUSSION

    Melampsora rusts on Salix are devastating pathogens for which the optimal control strategy is resistance breeding. An essential component to resistance breeding is the understanding of both host and pathogen biology. Although initial discoveries have been made in elucidating host resistance mechanisms (Cameron et al. 2008; Carlson and Gouker et al. 2019; Karp et al. 2011; Wilkerson et al. 2022), very little is known about the Melampsora rust population biology and diversity in North America. It is known that M. americana is the most abundant rust species on S. purpurea in this region (Bennett et al. 2011; Crowell et al. 2020; Kenaley et al. 2014; Smith et al. 2004) and that similar rust pathogens in Europe are generally panmictic (Barrès et al. 2008; Pei et al. 2007), yet we do not know the geographic distribution of any given M. americana genotype, nor do we have a sense of the diversity of these populations. Thus, the sampling scheme used here allows us to examine genetic diversity within a single shrub willow host across the Great Lakes region at specific sampling sites and across years at those sites.

    The lack of a reference genome for M. americana was a barrier to studies on pathogen biology and epidemiology. The generation of the genome sequence assembly of M. americana R15-033-03 version 1.0 provided a reference that enabled the discovery of SNP variants for population studies and will provide a resource for further research on willow leaf rust biology. Genome analyses revealed that M. americana has the greatest number of transporter proteins and transcription factors compared with the six other available annotated genomes from the genus Melampsora available on MycoCosm (https://mycocosm.jgi.doe.gov/). Additionally, the genome sequence of M. americana is the only one available for a willow rust pathogen. The addition of this genome sequence enables comparative genomics investigations, especially as they relate to host specificity differences between these plant pathogens.

    Our isolate collection strategies across the two years of the study were designed to address questions of geographic distribution and clonality. In 2015, samples were collected from a wide range of geographies, ranging from Vermont to Michigan and south to West Virginia, and factored in additional variables such as host domestication. Sample depth per location was low, and collections were opportunistic, because isolates could be collected only from diseased shrubs, leading to uneven sample depth between sites. In contrast, isolates collected in 2017 were limited to four cultivated S. purpurea fields in three locations. Sampling followed a more regimented scheme of 10 randomly distributed shrubs per field, and host genotypes were confined to variation observed within these breeding families. Additionally, each field was visited only once for a single collection timepoint. This sampling strategy enabled a broad overview of M. americana diversity in 2015, followed by a more in-depth analysis at a single field level in 2017. There was no sampling of M. americana in 2016 because of a lack rust observed in the field, probably because of hot, dry environmental conditions.

    To analyze the isolates collected in 2015, a geographic distance cutoff between any two isolates using hierarchical clustering of 80 km was used. Clustering resulted in 10 defined geographic groups, the majority of which fell within New York State. Neighbor joining analysis showed that there were no large monophyletic clades containing most isolates from a single geography. Instead, many smaller monophyletic clades containing isolates from one or multiple neighboring shrubs clustered together within large clades containing diverse isolate collection sites. This finding implies that immediate proximity of hosts may play a larger role in clustering of isolates than the a priori geographic groupings. Another interesting result is the evidence for large geographic spread of rust clonal lineages as defined by our GBS genotypes. An extreme example of this spread is an isolate collected from group 1 in Portland, New York clustered with an isolate in group 6 collected in Vermont (Fig. 2B), for a straight-line distance of approximately 480 km between collection sites. Such a distance between collection sites suggests that isolate gene flow may occur across large geographic distances. Such phenomena may play a role in the lack of population structure observed by principal component analysis, which revealed that there is no clustering of genotypes by geographic origin of collection. However, it is important to note that clonality was determined with a 96% identity cutoff, and the preceding example may be the result of type I error.

    The sampling scheme of isolates collected in 2017 enabled more directed questions about population differentiation by sampling site. Neighbor joining analyses of the four fields showed no large monophyletic clades containing all isolates from a single location. Instead, smaller clades containing isolates with matching field collection sites were part of larger monophyletic clades containing isolates from differing locations. This classification is supported by pairwise Fst and pairwise Jost D calculations for field sites Geneva, New York 1, Geneva, New York 2, and Portland, New York; however, the results did support moderate levels of population differentiation for the population in Morgantown, West Virginia. This matched the genetic diversity indices where the population in Morgantown, West Virginia displayed a lower level of genotypic diversity than the other three locations. This finding also aligns with the analysis of clonality, which identified one large clonal lineage included most isolates from Morgantown, West Virginia. The principal component analysis indicated that it is resolved from other genotypes by principal components 1 and 2. These data differentiate the Morgantown, West Virginia field from the other populations, yet there are genotypes from Morgantown, West Virginia persisting in multiple locations. DAPC was used to explore population structure with no a priori population designation, resulting in six identified clusters. Population structure with and without a priori population designation was compared in Supplementary Figure S4 and revealed highly similar patterns. This finding provides additional evidence of moderate population separation by geographic collection site.

    Taken together, both years of isolate collections support a hypothesis that M. americana has limited population separation by geography across the locations sampled. This hypothesis matches with what is known about a closely related rust species, M. larici-populina, on poplar in Europe (Barrès et al. 2008; Xhaard et al. 2011). Barrès et al. (2008) analyzed microsatellite and specific SNP marker genotypes of M. larici-populina isolates collected mostly from Europe to determine the Fst of geographically disparate populations. They found that there was virtually no distinction of populations across most of Europe, spanning >950 km. They did observe a slight shift in allele frequencies based on geographic distance, similar to what we saw in the isolates collected in Morgantown, West Virginia (group 9) in both 2015 and 2017. The only strong evidence for population separation was in the comparison of isolates from Canada and Iceland to those from Europe, for which large bodies of water separate potential hosts (Barrès et al. 2008). Further studies of M. americana in North America could include increased sampling across a broader geographic area to determine whether there are physical barriers that separate M. americana populations.

    In addition to the genetic diversity of pathogens within a field, the diversity of isolates infecting a single shrub or a single field is important for shrub willow breeding efforts. This information may yield insights about the expected success of deployed host resistance in a shrub willow cultivar in a field setting. To address this question, we attempted to collect up to three rust isolates per infected shrub in the field and compare frequency of clonality across these isolates. In many cases, only one isolate survived the process of culturing and uredospore collection. In total, we were able to isolate two or more single-spore isolates from 34 infected shrubs in 2015 and 29 shrubs in 2017 (Table 2). When two isolates were collected from a single shrub, a minority of them were clonal with (about 22% in 2015 and 25% in 2017). However, when three isolates were collected from the same shrub, about 42% were clonal in 2015 and 76% were clonal in 2017. The isolates collected in 2017 were collected exclusively from cultivated and densely planted willow fields, and it is possible that clonal spread via uredospores would be greater in this setting compared with individual wild shrubs with greater distances between hosts, as was the case in 2015. This information is useful because it may predict that disease within a single high-density field is likely to contain larger groups of clonal isolates. Although there may be greater influence of clonality within a field, it is important to note that all data collected in these experiments indicate that the pathogen population of this rust is highly diverse, and all hosts in the sampled regions are likely to be challenged by multiple genotypes of M. americana. This finding suggests that the best control methods may be to pyramid resistance genes in the hopes to provide durable resistance to this highly diverse pathogen population.

    Because M. americana is macrocyclic and heteroecious, it is expected that it sexually recombines between each shrub willow growing season on the alternate host, balsam fir. We sampled repeatedly in 2015 and 2017 for three sites. When we applied a 96% IBS clonality cutoff for all isolates together, it became clear that no genotype persisted across the two years of sampling in this study. However, this study was limited to two growing seasons and rarely repeated these samplings on the same host or field over these two nonconcurrent years. Expanding this study over many years, similar to Susi et al. (2020), who completed a large genetic study of M. lini over 16 annual epidemics in Australia, may provide greater evidence of the ability of M. americana to asexually overwinter.

    This investigation into the population of M. americana provides a foundation for future studies with broader geographic range and more intensive repetition of sampling. This study relied on sampling of infected Salix hosts once per growing season, at the end of the season in late summer to early fall. This method did not allow us to determine how the M. americana population changes over the course of the growing season. Future studies with repeated collections over time may determine whether the highly diverse population we observed in this study is similarly diverse over the growing season. Additionally, sampling from wider geographic ranges would yield a greater understanding of population separation across North America. Ultimately, these data support the hypothesis that M. americana in the Great Lakes region is highly diverse, with limited geographic population structure.

    ACKNOWLEDGMENTS

    We thank Rebecca Crust, Rebecca Wilk, Lauren Carlson, Dawn Fishback, Jane Petzoldt, Holly Lange, Garrett Giles, and Colin Day for outstanding technical assistance.

    The author(s) declare no conflict of interest.

    LITERATURE CITED

    Funding: Funding was generously provided by the U.S. Department of Agriculture National Institute of Food and Agriculture grants 2015-67009-23957, 2018-68005-27925, 2019-67011-29701, and 2019-67011-29698. Melampsora americana genome sequencing and assembly were accomplished under Joint Genome Institute Community Science Program CSP17 502972, supported by the Office of Science of the U.S. Department of Energy under contract number DE-AC02-05CH11231.

    The author(s) declare no conflict of interest.