Global Analysis of Hemileia vastatrix Populations Shows Clonal Reproduction for the Coffee Leaf Rust Pathogen Throughout Most of Its Range

Coffee (Coffea arabica and C. canephora) is one of the most important agricultural crops and one of the most popular beverages in the world (Giovannucci and Pierrot 2010; Spence and Carvalho 2020). It has an estimated retail market value of US$70 billion, and its production supports the livelihoods of millions of people in developing economies worldwide (Rhiney et al. 2021; Talhinhas et al. 2017). The earliest records of human consumption of coffee date to Ethiopia ca. AD 757 to 850 (Amamo 2014; Herrera and Lambot 2017; Rodrigues et al. 1975). Coffee was first commercially cultivated in Yemen and traded by merchants through much of the Arabian Peninsula in the 15th century (Herrera and Lambot 2017; Hindorf and Omondi 2011). In response to increased European and, later, North American demand for the beverage, vast coffee plantations were established in Southeast Asia in the 1600s, and by the 1800s coffee was being produced in much of the Americas (Herrera and Lambot 2017). Today, outside the major production centers of Brazil, the crop is grown primarily by small-shareholder farmers throughout the global coffee belt (between the tropics of Capricorn and Cancer) (Voora et al. 2019). Throughout its range, coffee diseases are a major limiting factor to production, the most important of which is coffee leaf rust (CLR), caused by the fungus Hemileia vastatrix Berk. & Broome (Berkeley and Broome 1869; Hindorf and Omondi 2011; McCook 2006; Waller et al. 2007).

The first recorded outbreaks of CLR date to India and Sri Lanka (then Ceylon) in the late 1800s, where H. vastatrix was responsible for the collapse of the coffee plantation industry in Ceylon within 10 years of its arrival (McCook 2006; McCook and Vandermeer 2015). With the discovery of H. vastatrix in Hawaii in late 2020, CLR has now colonized all major coffee growing regions of the world (EPPO Global Database; McCook 2006; State of Hawaii Dept. of Agriculture 2020), causing reduced yield and marketability and threatening the livelihoods of subsistence farmers (Belachew et al. 2015; Hindorf and Omondi 2011; Nelson 2008; Rhiney et al. 2021). Today, global losses to CLR are estimated at US$1 to 2 billion per year (McCook 2006).

Breeding for rust resistance is the main control strategy, but gradual breakdowns of resistance have been documented in many improved coffee varieties in several countries (Reuben and Mtenga 2012; Rodrigues et al. 1975; Zambolim 2016). Recent resurgence of epidemic levels of disease in parts of the Americas between 2008 and 2013 have renewed interest in developing effective control measures for the fungus (Avelino et al. 2015; Cristancho et al. 2012). Resistance breakdown could be caused by changes in pathogen populations and emergence of new physiological races (Kiyosawa 1982; McDonald and Linde 2002). The emergence of new physiological races of CLR, in turn, may arise by different means that include the process of natural selection acting on populations that have been exposed to increased deployment of individual resistant cultivars or the result of new combinations of avirulence genes through sexual reproduction (Burdon et al. 2016; Kiyosawa 1982). Because the mechanisms behind the emergence of new races depends on the reproductive strategy of the pathogen (Kiyosawa 1982), understanding pathogen biology is essential to effective deployment of new cultivars and management strategies.

H. vastatrix is an obligate biotroph belonging to the order known collectively as rust fungi (Pucciniales, Basidiomycota) (Aime et al. 2006). Many rust fungi alternate between two unrelated host plants in order to complete their life cycle, called heteroecism (Aime et al. 2018; Littlefield 1981; Voegele et al. 2009). One host is colonized by the fungal gametothallus and is the site of sexual recombination via fertilization (Aime et al. 2018; Kolmer et al. 2009). The fungal sporothallus colonizes an unrelated host, where the asexual repetitive portion of the life cycle occurs via continued clonal propagation of urediniospores (Aime et al. 2018). Examples of heteroecious rust include the wheat rust Puccinia graminis, which alternates between hosts in the Berberidaceae (gametothallus) and Poaceae (sporothallus) (Leonard and Szabo 2005).

Some rust fungi, called autoecious, have adapted to complete both parts of the life cycle on a single host that is colonized by both the sporothallus and the gametothallus. Examples include the sunflower rust, Puccinia helianthi (Sendall et al. 2006) and the pathogen of Russian knapweed, P. acroptili (Bruckart et al. 2010). Because of the complexity of rust life cycles and extreme morphological variation between gametothalli and sporothalli, complete life cycles are unknown for thousands of rust fungal species, including some of agricultural importance (Aime and McTaggart 2021). For example, the complete life cycle of the wheat rust pathogen P. striiformis was not discovered until 2010 (Jin et al. 2010; Zhao et al. 2013), despite >150 years of study on this pathogen. In the case of H. vastatrix only the sporothallus is known, which infects several Coffea species (Coutinho et al. 1995; Talhinhas et al. 2017). Establishment of a gametothallus and fertilization have never been observed on Coffea spp. or on any other host species (Aime and McTaggart 2021; Coutinho et al. 1995). It remains unknown whether a gametothallus is formed by this rust and, if so, whether it is formed on an unknown alternate host or on the same host as the sporothallus.

Deciphering the extent and capacity of a fungal pathogen to undergo sexual recombination can have enormous consequences for disease management and breeding programs. For example, sexually reproducing species are often able to overcome a host’s resistance more efficiently than those with only asexual reproduction (McDonald and Linde 2002). One example of this difference is the new wheat stem rust (P. graminis f. sp. tritici) race, Ug99, which is believed to have arisen in central-eastern Africa (Pretorius et al. 2000). The new combination of genes that came together in Ug99 has been capable of overcoming resistance in almost every developed wheat cultivar (Aydoğdu and Boyraz 2012; Singh et al. 2011). And, as previously mentioned, more virulent races in clonally propagating fungi can emerge under intense deployment of resistant host varieties, which can lead to an erosion of host resistance through directional selection (Burdon et al. 2016).

Several indirect approaches involving analyses of population-level variation use molecular markers to infer the reproductive strategies of plant pathogens (Nieuwenhuis and James 2016; Taylor et al. 1999). Such methods remain the best means for inferring reproduction in obligate pathogens, such as rust fungi, that cannot be experimentally manipulated in culture. Various molecular markers have been applied to the question of reproductive mode in H. vastatrix, such as amplified fragment length polymorphisms (AFLPs) (Cabral et al. 2016), random amplified polymorphic DNA (RAPDs) (Cabral et al. 2016; Gouveia et al. 2005), and single-nucleotide polymorphisms (Silva et al. 2018), with conflicting results. The use of one type of neutral marker, simple sequence repeat (SSR) markers, is useful in population genetics because of their neutrality, repeatability, and high variability that make them superior to many other types of markers (Kraic et al. 1998; Powell et al. 1996; Rawat et al. 2014; Vieira et al. 2016). Additionally, because SSRs need low quantities of DNA template, these markers are ideal for analyses of older herbarium specimens or obligate pathogens that are typically embedded in host tissues (Schuelke 2000; Senan et al. 2014). Analysis of SSR data has been used to resolve population structure and outcrossing potential for several rust fungi such as P. graminis f. sp. tritici, P. striiformis, P. striiformis f. sp. tritici, P. triticina, P. psidii, and Melampsora larici-populina (Ali et al. 2014; Cheng and Chen 2014; Goyeau et al. 2007; Graça et al. 2013; Pernaci et al. 2014; Thach et al. 2016; Visser et al. 2019).

The goal of the present study was to analyze the global population structure of H. vastatrix to infer whether the fungus is undergoing random mating, which would indicate sexual recombination. To do this, we developed a set of eight SSRs from the H. vastatrix genome and assembled a global set of 105 H. vastatrix isolates collected from coffee plantations in both hemispheres between 2014 and 2019. H. vastatrix isolates were grouped based on clustering analysis, and a combination of statistical approaches was applied to infer the reproductive mode within studied H. vastatrix populations.

MATERIALS AND METHODS

H. vastatrix specimens.

A total of 105 isolates of H. vastatrix were collected from around the world between 2014 and 2019 (Supplementary Table S1). All specimens were collected on C. arabica from commercial coffee plantations and farms except for three specimens collected from naturalized plants in the Darien Province of Panama (MCA8132_HV1-3) and the specimens from Cameroon, which were collected on C. canephora plants from a small village where coffee is grown in householder gardens for local consumption. Each collection consisted of infected leaves taken from a single host tree and accompanied with metadata on host variety, locality, and pertinent ecological conditions. Coffee leaves were dried in the field in a plant press and deposited at the Arthur Fungarium at Purdue University.

Genome assembly and SSR primer design.

Approximately 1.5 µg of genomic DNA was extracted from H. vastatrix MCA6123 (Portland Parish, Jamaica, 26 January 2016) from urediniospores with a PowerPlant Pro DNA isolation Kit (MoBio, Carlsbad, CA) according to the protocol of Aime (2006). Library preparation for whole genome sequencing was done with Illumina’s TruSeq DNA PCR-Free Prep Kit (Illumina, San Diego, CA) at the Purdue Genomics Core Facility (PGCF). Sequencing was performed at the PGCF on one quarter lane of an Illumina HiSeq2500 platform. The genome assembly of H. vastatrix was performed via a de novo approach that assumes no prior knowledge of the source of sequence. The paired-end reads were assembled with ABySS version 1.9.0 with a k-mer length of 70 (Simpson et al. 2009). Genome assembly statistics were generated in QUAST version 5 (Mikheenko et al. 2018). Our assembly was used as the input for the microsatellite detection and primer design with QDD version 3.2.1 (Meglécz et al. 2009, 2014), according to the approach of Díaz-Valderrama and Aime (2016). The selection criteria for primer design were as follows: number of motif repeats of 5 to 30, product size of approximately 186 bp, primer length of approximately 21 nucleotides, and primer melting temperature of 57 to 60°C.

DNA extraction.

For each isolate (Supplementary Table S1), a single sorus (about 3 × 3 mm) from a single leaf was extracted according to standard protocols with a DNeasy PowerPlant Pro Kit (Qiagen, Hilden, Germany) or PowerPlant Pro DNA Isolation Kit (MoBio, Carlsbad, CA). To avoid specimen crossover, sterile razor blades and gloves were changed and the work surface sterilized with ethanol or DNA Away (Thermo Fisher Scientific, Waltham, MA) between excisions. In a few cases where the sori were too small to obtain sufficient DNA, two or three adjoining sori from a single leaf and infection ring were excised. DNA was quantified with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific) and diluted to 0.2 ng/μl before SSR amplification. Identity of extracted DNA was confirmed by sequencing of the nuclear large subunit rDNA according to protocols of Aime (2006) for specific isolates showing unusual genotype patterns to confirm identity as H. vastatrix.

Microsatellite marker screening and reproducibility.

Screening of potential microsatellite loci was performed via a modified M13-tailed primer method (Schuelke 2000), as previously described (Díaz-Valderrama and Aime 2016; Koch and Aime 2018). Initially we screened a previously published genome (Cristancho et al. 2014) for potential microsatellite regions but were unsuccessful.

From our newly generated genome we identified 121 microsatellite markers located on separate scaffolds (Supplementary Table S2) from our genome assembly. These were initially screened against a panel of six H. vastatrix specimens collected from Puerto Rico, Jamaica, and Cameroon. Fragment amplification, sequencing, and analysis followed methods outlined in the next section. Of the initially screened 121 SSR markers, 113 were monomorphic across all samples. The remaining eight showed polymorphic bands and were selected for full genotyping (Supplementary Table S3).

We tested the reproducibility for all SSR markers by randomly selecting four H. vastatrix isolates and generating three full sets of genotypes for each isolate at an annealing temperature of 49°C. With each isolate displaying identical genotypes in each of the three repetitions, we proceeded with genotyping of the entire dataset. To ensure consistency, each round of amplifications also included one specimen of H. vastatrix with a known genotype as a positive control, as well as a negative control consisting of reagents plus sterile water instead of DNA template.

SSR amplification and fragment analysis.

For the microsatellite loci, reactions were performed at a volume of 12.5 μl, with 6.25 μl of Promega 2× Taq DNA Polymerase PCR Master Mix (Promega Co., Madison, WI), 0.16 μl of the forward primer with the M13 tail, 0.47 μl of the M13 primer with one of four fluorescent dyes (6-FAM, NED, PET, and VIC, Thermo Fisher Scientific, Waltham, MA), 0.63 μl of the reverse primer (all at a concentration of 10 μM), and 5 μl of DNA template. The PCR conditions standardized for microsatellite amplification in this study were 94°C for 5 min; 30 cycles of 94°C for 30 s, 49°C for 45 s, and 72°C for 45 s; eight cycles of 94°C for 30 s, 53°C for 45 s, and 72°C for 45 s; and 72°C for 10 min. Amplification of the eight markers was confirmed by running the amplicons on 2% agarose gel in 1× TAE electrophoresis buffer, stained with Gel Red, visualized with a Gel Doc EZ Imager (Bio-Rad, Hercules, CA), and printed with a Digital Graphic Printer UP-D897 (Sony Corporation of America). All gels were run at 80 V constant power supply for 45 min. After confirmation, PCR products were diluted in nuclease-free water (1:100), of which 2 μl was added into 10 μl containing 1 ml of HiDi formamide loading buffer (Thermo Fisher Scientific) and 12.5 μl of the GeneScan 500 LIZ dye size standard (Thermo Fisher Scientific). Diluted PCR products were sent for fragment analysis at the PGCF with an ABI 3730XL sequencer used to generate fragment data. The fragment sizes for each sample at each locus were analyzed in Geneious version 9.1.8 (Biomatters Ltd., Auckland, New Zealand) (Kearse et al. 2012).

Population and random mating analyses.

A global data matrix was constructed from all isolates (n = 105) for which complete data at all eight SSR markers were obtained (Supplementary Table S4). The data matrix was transformed to GenAlEx version 6.5 format (Supplementary Table S4) (Peakall and Smouse 2012, 2006). All analyses were performed with R package Poppr version 2.8.3, unless otherwise specified. To determine the minimum number of loci necessary to discriminate between genotypes, we generated a genotype accumulation curve, randomly sampled 1,000 times for each box plot (Kamvar et al. 2015). Multilocus genotypes (MLGs) were computed to determine the evenness of genotype distribution between populations (Kamvar et al. 2014).

Because the statistical power to infer reproductive mode requires large sample sizes (Brown 1975; Grünwald et al. 2017; Milgroom 1996), we initially grouped our H. vastatrix isolates into three large predefined populations for these analyses: global (all specimens, n = 105); New World (Bolivia, Colombia, El Salvador, Guatemala, Honduras, Jamaica, Panama, Peru, and Puerto Rico, n = 62); and Old World (Cameroon, Congo, and India, n = 43). We also analyzed additional populations that resulted from the Structure analysis: cluster 1, called Robusta; cluster 2, called Arabica; and cluster 3, called Robusta/Arabica. All datasets were complete, and all analyses were performed with and without clone correction (Grünwald et al. 2017; Kamvar et al. 2017, 2014).

To explore the population structure of H. vastatrix and determine the most likely number of clusters to which our isolates belong, we used the method of assignment implemented in Structure (Pritchard et al. 2000). About 30 independent Markov chain Monte Carlo replicates of 100,000 steps were run after a burn-in period of 20,000 steps (Ali et al. 2014; Jombart et al. 2010). The structure outputs were processed with Structure Harvester version 0.6.94 (Earl and vonHoldt 2012) to determine the optimal K and CLUMPAK version 1.1.2b (Kopelman et al. 2015) to determine cluster membership coefficients for each isolate, then visualized on Structure Plot version 2.0 (Ramasamy et al. 2014). Population differentiation was also estimated between clusters via Hedrick’s G_ST estimator implemented in Modern Methods of Differentiation (mmod) package version 1.3.3 (Hedrick 2005). It is presumed that values close to 0 do not differentiate between H. vastatrix populations, whereas values close to 1 indicate differentiation between populations.

The population allelic richness (AR) and private allelic richness (pAR) for the eight tested markers were computed by the rarefaction method implemented in HP-RARE 1.0 software (Kalinowski 2005). To infer a reproductive mode for H. vastatrix, we performed four analyses as follows. First, we used the index of association (I_A) and the standardized index of association (rbarD) among H. vastatrix populations (Grünwald et al. 2017; Kamvar et al. 2017, 2014). The assumption of random mating (where linkage between loci is not expected) is rejected when the observed value for rbarD falls outside the bell curve produced by the simulation of a randomly mating population (Agapow and Burt 2001). For this analysis, we used 1,000 random permutations for all populations we tested. A pairwise rbarD over the eight loci was used to ensure that the pattern of linkage disequilibrium observed was not caused by signal from a single pair of loci.

Second, Hardy–Weinberg equilibrium (HWE) of H. vastatrix populations was analyzed in GenAlEx version 6.5 (Peakall and Smouse 2012, 2006). HWE was computed via the χ² test over H. vastatrix populations based on the expected genotype frequencies calculated from the allelic frequencies. For our eight loci we tested the hypothesis that genotype frequencies do not deviate from HWE. Populations under HWE are expected to be undergoing sexual recombination; absence of HWE is expected in clonally reproducing populations. Basic genetic parameters including number of different alleles (Na), effective number of alleles (Ne), and observed (Ho) and expected (He) heterozygosity were calculated across loci for each population. The value of these measures will range from zero (no heterozygosity) to nearly 1.0 (for a system with a large number of equally frequent alleles).

Third, clonal richness [(the number of genotype − 1)/(number of isolates − 1)] was calculated for each H. vastatrix population (Arnaud-Haond et al. 2007). It is assumed that values close to 0 are typical of a strictly clonal population, in which all individuals share the same genotype, whereas a value close to 1 suggests that all individuals in the population have distinct genotypes, as expected under sexual reproduction.

Finally, we used the interactive tool imsn() with 1,000 random seeds to generate a minimum spanning network (MSN) based on the eight SSR markers. An MSN allows visualization of genetic relatedness between individual MLGs in the global H. vastatrix population (Kamvar et al. 2014, 2015). Populations with sexual recombination will show patterns of reticulation in the MSN, in contrast to populations undergoing clonal reproduction (Linder et al. 2003). We used Nei’s genetic distance as it has been applied in other analyses of rust fungal populations (Olivera et al. 2015).

RESULTS

SSR markers.

The resulting 333,481,311-bp draft assembly of the H. vastatrix genome from Jamaica comprises 72,027 contigs >1,000 bp and has an N₅₀ of 2,601 bp and an L₅₀ of 24,839 contigs. Within this H. vastatrix genome, we screened 121 SSR markers from different scaffolds that were consistently amplifiable across an initial tester set of six H. vastatrix isolates collected from three geographic regions (Puerto Rico, Jamaica, and Cameroon). Of these, only eight were found to be polymorphic and retained for additional genotyping of all 105 H. vastatrix isolates. Each of the eight selected SSR markers was found to be reproducible at an annealing temperature of 49°C (Supplementary Table S3). The genotype accumulation curve shows that across levels loci analyzed were sufficient to observe the total number of MLGs, although the addition of extra markers would probably increase the number of recovered genotypes in H. vastatrix populations (Supplementary Fig. S1).

MLGs.

A total of 12 MLGs were found (MLG 1 to MLG 12), distributed among the 105 H. vastatrix isolates (Fig. 1 and Supplementary Fig. S2). We found unique MLGs that are present in only one country: Cameroon (MLG 12; n = 4), Congo (MLG 3; n = 1), Guatemala (MLG 5; n = 1), El Salvador (MLG 6, MLG 11; n = 1 for each MLG), and India (MLG 1 and MLG 4 [n = 1 for each MLG], MLG 7 and MLG 8 [n = 5 and 6 for each MLG, respectively]) (Fig. 1 and Supplementary Fig. S2). Some MLGs were shared between some but not all countries: Congo, Colombia, and Puerto Rico (MLG 2; n = 31); India and Colombia (MLG 9; n = 2 for each country); and Bolivia, El Salvador, Guatemala, Honduras, India, Jamaica, Panama, Peru, and Puerto Rico (MLG 10; n = 49) (Fig. 1 and Supplementary Fig. S2). One MLG was represented by only four isolates belonging to Cameroon (MLG 12), and 11 MLGs were shared between the remaining countries (Fig. 1 and Supplementary Fig. S2).

Download as PowerPoint

Tests for genetic structure.

Our data partitioned the global population of H. vastatrix into three clusters (Fig. 2): cluster 1, herein called Robusta, consisted solely of the four Cameroonian isolates, all of which were collected on C. canephora grown for local consumption; cluster 2, herein called Arabica, contains specimens from Congo (n = 23), Colombia (n = 8), India (n = 1), and Puerto Rico (n = 1); and cluster 3, herein called Robusta/Arabica, contains all other analyzed isolates (Colombia, n = 2; Bolivia, n = 1; El Salvador, n = 8; Guatemala, n = 9; Honduras, n = 2; India, n = 15; Jamaica, n = 9; Panama, n = 13; Peru, n = 8; Puerto Rico, n = 1). However, Hedrick’s G_ST values indicate a very low level of genetic differentiation between Arabica and Arabica/Robusta (non-clone-corrected dataset, pairwise G_ST = 0.204; clone-corrected dataset, pairwise G_ST = –0.038), with only Robusta showing genetic differentiation from the other two clusters (non-clone-corrected dataset, pairwise G_ST = 0.841 and 0.786; clone-corrected dataset, pairwise G_ST = 0.631 and 0.625).

Download as PowerPoint

The AR was similar for the global (AR = 3.5), Old World (AR = 3.38), Robusta/Arabica (AR = 2.5), and New World (AR = 2.29) populations of H. vastatrix but lower in the Arabica (AR = 2) and Robusta (AR = 1.75) populations (Table 1). The pAR, determined by the rarefaction method, ranged from 0.09 (New World) to 3.5 (global) (Table 1). Similar results were obtained from the clone-corrected dataset for all populations (Table 1).

TABLE 1. Population analyses in Hemileia vastatrix^a

Population	N	MLG	IA	p.IA	rbarD	p.rbarD	Clonal richness (MLG – 1)/(N − 1)	Average of allele richness
								AR	pAR
Non-clone-corrected
Global	105	12	2.735	0.001	0.434	0.001	0.106	3.5	3.5
New World	62	6	0.604	0.001	0.241	0.001	0.082	2.29	0.09
Old World	43	9	3.602	0.001	0.557	0.001	0.19	3.38	1.17
Robusta	33	3	NA	NA	NA	NA	0	1.75	0.88
Arabica	101	11	−0.031	0.613	−0.031	0.989	0.1	2	2
Robusta/Arabica	68	8	−0.067	0.904	−0.025	0.999	0.104	2.5	2.5
Clone-corrected
Global	12	12	2.303	0.001	0.353	0.001	1	3.5	3.5
New World	6	6	−0.508	0.926	−0.173	1	1	2.29	0.09
Old World	9	9	2.661	0.001	0.409	0.001	1	3.38	1.17
Robusta	1	1	NA	NA	NA	NA	NA	1.75	1.12
Arabica	3	3	−0.500	0.661	−0.500	0.999	1	2	2
Robusta/Arabica	8	8	−0.592	0.998	−0.198	1.000	1	2.5	2.5

^a AR, allelic richness; I_A, index of association; MLG, number of multilocus genotypes observed; N, number of individuals sampled; NA, not applicable; pAR, private allelic richness; p.I_A, P values of the index of association; p.rbarD, P values of the standardized index of association; rbarD, standardized index of association.

View as image HTML

Tests for sexual reproduction.

A hypothesis of random mating (i.e., sexual outcrossing) in the global population of H. vastatrix was rejected in global analyses of both clone-corrected (I_A = 2.303, P = 0.001; rbarD = 0.353, P = 0.001) and non–clone-corrected (I_A = 2.735, P = 0.001; rbarD = 0.434, P = 0.001) datasets (Table 1 and Supplementary Fig. S3); within populations defined as Old World and New World (Table 1); and within clusters detected by Structure: Arabica (clone-corrected data: I_A = −0.500, P = 0.661; rbarD = −0.500, P = 0.999; non-clone-corrected data: I_A = −0.031, P = 0.613; rbarD = −0.031, P = 0.989), Robusta/Arabica (clone-corrected data: I_A = −0.592, P = 0.998; rbarD = −0.198, P = 1; non-clone-corrected data: I_A = −0.067, P = 0.904; rbarD = −0.025, P = 0.999). Robusta was not tested because of its small sample size (n = 4) (Table 1).

Of the eight SSR loci tested, six deviated from HWE (P ≤ 0.05), and the other two (154 and 180) were inconclusive (Supplementary Table S5). In clone-corrected datasets the opposite was found, where some loci do not deviate from HWE (P ≥ 0.05). The percentage of polymorphic loci ranged from 75 to 100% for all analyzed populations. The mean Ho ranged from 0.75 to 0.87; He ranged from 0.38 to 0.51. In clone-corrected datasets these values were Ho = 0.75 to 0.86 and He = 0.38 to 0.55, respectively (Supplementary Table S5). The highest proportion of observed and expected heterozygosity was detected in the global, New World, Old World, and Robusta/Arabica populations (global and New World Ho = 0.83, He = 0.48; Old World Ho = 0.80, He = 0.51; Robusta/Arabica Ho = 0.87, He = 0.45) (Supplementary Table S5). The lowest mean observed and expected heterozygosity ratio was found in the Structure-defined populations (Robusta Ho = 0.75, He = 0.38; Arabica Ho = 0.75, He = 0.38) (Supplementary Table S5). The same pattern was showed in the clone-corrected dataset (Supplementary Table S5). Clonal richness is low in all datasets (global = 0.106, New World = 0.082, Old World = 0.19, Robusta = 0, Arabica = 0.1, and Robusta/Arabica = 0.104) (Table 1). All populations except Robusta have a clonal richness of 1 for the clone-corrected datasets (Table 1). Finally, MSN analysis lacks reticulation in the terminal nodes, suggestive of a clonal pattern of reproduction for the global H. vastatrix population (Fig. 1).

DISCUSSION

Herein we developed a set of SSR markers for H. vastatrix and used them to screen a global sample of 105 H. vastatrix samples isolated primarily from C. arabica coffee plantations across 12 coffee producing countries. Clonally reproducing fungi have lower genetic diversity than sexually reproducing fungal populations and are thereby less enriched in polymorphic SSR markers (Grünwald et al. 2017). In stark contrast to the comparatively rich percentage of informative markers developed from primarily sexually reproducing fungi (Koch and Aime 2018) and even compared with other clonally propagating plant pathogens (Díaz-Valderrama and Aime 2016), we were able to identify only eight polymorphic markers of 121 screened, which may be the first indicator of a primarily clonal strategy for this pathogen. Selected markers were derived from different scaffolds within the H. vastatrix assembled genome, reducing the likelihood of any linkage between them.

A total of 12 MLGs (12 MLG/n = 105, or 0.11) were identified from this sampling, which is similar to the ratio found for a clonal population of the rust fungus P. striiformis f. sp. tritici in the United States (MLGs = 32/n = 270, or 0.12) (Cheng and Chen 2014). However, the genotype accumulation curve demonstrates that we were able to discriminate between all MLGs, within a margin of error, with the markers we applied (Supplementary Fig. S1). The distribution of MLGs in the H. vastatrix populations shows sharing of some MLGs, such as MLG 10, between almost all sampled countries (Fig. 1 and Supplementary Fig. S2). The existence of identical MLGs in different countries and regions also suggests that the global population of H. vastatrix is predominantly clonal (Anderson and Kohn 1998). The existence of identical genotypes (MLGs) in different geographic locations has been attributed to founder effects during historic or recent epidemic events in plant pathogens of monoculture crops where genotypes are rapidly propagated without losing effective combinations of genes (Drenth et al. 2019; McTaggart et al. 2020), as is the case for the urediniospore stage of rust pathogens (Aime et al. 2018). Successful genotypes in coffee plantations can be further spread via the bridgehead effect, wherein particularly successful invasive genotypes are the source of colonization into remote new territories (Lombaert et al. 2010). Either of these effects could be facilitated by effective long‐distance dispersal of urediniospores through strong winds or by infected host germplasm exchange (Talhinhas et al. 2017). The distribution of MLGs shows several unique MLGs in both hemispheres. Such a pattern may arise from local mutation events, following invasion into a naive population, as has been suggested for one subpopulation of H. vastatrix in an analyzed Brazilian population (Cabral et al. 2016). Additionally, MLG 10 consists of rust infections from both C. arabica and C. canephora. A similar result was previously reported from both Coffea species and derivatives in Híbrido de Timor and Icatu, in the major producing states of Brazil (Cabral et al. 2016). These studies suggest that successful H. vastatrix genotypes have the capacity to infect more than one Coffea species.

The MSN analysis (Fig. 1) also indicates a pattern of clonal reproduction and the bridgehead effect, with little or no reticulation. The reticulation that is present is caused by the uncertainty about the placement of the Cameroonian isolates, because most isolates belonging to only two MLGs from which smaller MLGs are derived (Linder et al. 2003; Lombaert et al. 2010).

Our data were treated and analyzed in different datasets: a global population that included all samples; regional populations consisting of an Old World dataset (Congo, Cameroon, and India) representing the original spread of H. vastatrix, which began in about 1860; a New World dataset (Latin American and Caribbean isolates) into which the pathogen began spreading in 1970 (McCook 2006); and three cluster-detected populations (Robusta, Arabica, and Robusta/Arabica) as determined by Structure analyses (Fig. 2) and differentiated via Hedrick’s G_ST. Our results reject a null hypothesis of random mating in the sampled H. vastatrix populations for all non-clone-corrected datasets, and six of eight markers tested in H. vastatrix populations deviate from HWE (P ≤ 0.05; Supplementary Table S5). This is in marked contrast to patterns found for the sexually reproducing rust P. striiformis f. sp. tritici, where Asian (Himalayan) populations of the wheat stripe rust did not depart from HWE, suggesting the occurrence of random mating within those populations (Ali et al. 2014).

Ho was consistently higher than He (Supplementary Table S5), meaning that any departures from HWE could be associated with an isolation-breaking effect rather than random mating (Hartl and Clark 1997). Significant (P ≤ 0.05, Supplementary Table S5) differences between the averages of He and Ho heterozygosity for all populations indicate clonality. Moreover, the low levels of genetic diversity (He = 0.38 to 0.51; Supplementary Table S5) and the detection of two monomorphic loci (154 and 180) because of the presence of one allele within some populations (New World, Robusta, Arabica, and Robusta/Arabica) indicate low genetic diversity within each analyzed H. vastatrix population, also supporting a model of clonal reproduction for this fungus. Considering the broad geographic sampling in our study, the amount of genetic diversity we found within sampled isolates of H. vastatrix is much lower than expected but still higher than that reported in previous studies with different markers and restricted geographic scope (Cabral et al. 2016; Gouveia et al. 2005; Nunes et al. 2009).

The average AR and pAR are highest in the Old World population, which is consistent with the known center of origin of the pathogen residing in the Old World (Mboup et al. 2009; Saleh et al. 2014). The low clonal richness in all datasets indicates that random mating has not played a role in the extant sampled H. vastatrix populations. These values are similar to those found for the wheat rust P. striiformis f. sp. tritici in the United States, where the clonal richness was too low (0.14 and 0.17) to be considered a sexually reproducing population (Cheng and Chen 2014). Recently, the Himalayan and near-Himalayan regions have been pinpointed as the center of diversity for P. striiformis f. sp. tritici, and this region has been shown to harbor populations that do not depart from HWE (Ali et al. 2014), in contrast to what was found in North America for that rust. The MSN analyses of our H. vastatrix samples (Fig. 1) are also inconsistent with a pattern of random mating, as is seen in sexually reproducing populations such as the blueberry fungus pathogen Exobasidium maculosum (Stewart et al. 2015) or the rust fungus Austropuccinia psidii (McTaggart et al. 2020).

Our data indicate that most sampled H. vastatrix isolates are derived from just two clones: MLG 2, which is found in the Congo, Colombia, and Puerto Rico, and MLG 10, which was found in India and in all New World countries except Colombia (Fig. 1 and Supplementary Fig. S2). MLG 10 most closely represents the known historical spread of the disease from Asia into the Americas; whether MLG 2 represents a separate invasion into the New World from Africa will require additional sampling. Of the countries sampled, India contains the greatest number of MLGs (Supplementary Fig. S2), consistent with the fact that India and Sri Lanka saw the first outbreaks of epidemic H. vastatrix in the 1870s (McCook 2006) and from which MLG 10 probably spread because of the founder and bridgehead effects already discussed. The sampled population in Cameroon is distinct from the rest of the sampled global population (Figs. 1 and 2, and Supplementary Fig. S2). Whether this difference reflects the low sample size from this country (n = 4) or is long isolation and selection on a different host (C. canephora) remains to be determined.

Although our study is the first to analyze a statistically significant global population of H. vastatrix, our sampling is still incomplete in two notable regions, Brazil (the single largest producer of coffee) and north central Africa (the origin of cultivated coffee). Brazil is the first documented country in the New World to be invaded by H. vastatrix (McCook 2006), and sexual recombination has been inferred in a subset of two of five tested populations of H. vastatrix in Brazil through AFLPs (Cabral et al. 2016). Although we were unable to detect a signal of sexual reproduction in the New World, we cannot rule out the possibility that the Brazilian population acts as a center for diversity from which clonal lineages have been dispersed elsewhere into the Americas. However, the low sampling of H. vastatrix in the cited study (n = 10 from Paraná and n = 20 from São Paulo; Cabral et al. 2016) could also result in erroneous interpretations. For example, in our study after sampling was grouped by region, we found we could not reject a hypothesis of sexual recombination in South America alone because of the small sample size (n = 19) (non–clone-corrected data: I_A = 0, P = 0.373; rbarD = 0, P = 0.379). In further analyses of the South American dataset, we found that the pairwise rbarD over the tested eight loci showed that the signal of linkage disequilibrium was caused by various pairs of monomorphic loci (SSR148, SSR150, SSR159, SSR166, SSR180) containing undefined values and not a true signal of recombination. Of the Brazilian populations defined in Cabral et al. (2016), the small sampling of São Paulo and Paraná and the lack of pairwise rbarD analyses could influence interpretation of reproduction status in both regions; when the dataset is analyzed with the larger sample size (n = 112), reproduction was found to be clonal Cabral et al. (2016).

The other major region missing from our analysis is north central Africa, where wild coffee plants, and presumably H. vastatrix, are believed to have originated (Herrera and Lambot 2017). In a previous study that applied single-nucleotide polymorphisms on Old World populations that included multiple isolates from India, Tanzania, Kenya, Philippines, and Timor and a single isolate each from Brazil and Ethiopia, a signal for recombination was detected in one domesticated lineage of H. vastatrix (Silva et al. 2018). Previous studies have shown that the center of origin of rust pathogens is generally associated with the center of diversity, which often shows strong evidence for sexual recombination (Mboup et al. 2009; Saleh et al. 2014). And there is precedent for rust fungi, such as the examples outlined previously for P. striiformis f.sp. tritici (Ali et al. 2014; Cheng and Chen 2014), that reproduce sexually near the point of origin, presumably because of the availability of codistributed gametothallus hosts, but reproduce solely as clones outside this range when adequate sporothalli hosts exist, such as in cultivated crops, to ensure constant cycling in absence of the gametothallus host (Aime et al. 2018). However, most population geneticists recommend analyses of data both before and after clone correction to provide a better understanding of the influences of each reproductive mode (clonal and sexual reproduction) on population structure (Grünwald et al. 2017), and whether the data were clone corrected or not was not specified in the study by Silva et al. (2018). Although we were unable to detect a signal of sexual reproduction in the Old World, this population does contain a greater number of MLGs (9 of the 12 MLGs vs. 6 of 12 for the New World population) that indicates populations here are more diverse than those dispersed to the New World.

Other studies with sexually reproducing rust fungi such as P. striiformis f. sp. tritici (n = 19 to 20 SSRs) and P. graminis f. sp. tritici (n = 17 SSRs) (Ali et al. 2014; Thach et al. 2016; Visser et al. 2019) were able to generate more SSR markers than were used in the present study, and we cannot rule out the possibility that eight loci could be insufficient to detect a signal of sexual recombination in H. vastatrix (Grünwald et al. 2003, 2017). Further analysis increasing the number of samples per location and using a larger set of SSR polymorphic markers (e.g., McVean et al. 2004) in H. vastatrix would improve the strength of the inferences drawn from our data.

The global development of CLR disease caused by H. vastatrix has been a challenge for farmers and researchers since the first epidemic outbreak in 1869 in Ceylon and southern India (McCook 2006). Using multiple lines of inference, we were unable to detect any signal of random mating for this fungus within a global dataset of H. vastatrix isolates, suggesting that the fungus reproduces clonally throughout its invaded range, wherever sufficient populations of host plants exist to support its growth. Our findings suggest that the process of directional selection, rather than recombinant genotypes, may be an important driver of observed erosion of host resistance.

ACKNOWLEDGMENTS

Collecting permits were issued by Cameroon’s Ministry of Research and Scientific Innovation; Ing. Juan José Aparicio Porres, La Paz, Bolivia; the Peruvian National Forestry and Wildlife Authority (SERFOR); and collection permit SE/PH-3-16 and associated export permit from the Ministry of the Environment of Panama (MiAMBIENTE). Additional samples were collected by the Inter-American Institute for Cooperation in Agriculture and World Coffee Research in partner countries (E. Johnson, S. Pruvot-Woehl, and C. Montagnon). We thank Phillip SanMiguel at the Purdue University Genomics Core for assistance with genome sequencing. Finally, we are grateful to Purdue University undergraduates Cheyanne Marie Woolwine, Thijs Montalvo, and Matthew Alexander Schwarzkopf for assistance with laboratory work.